APL/JHU 
TG  1293 A 
JULY  1976 
Copy  No. 


00  • 
00 
05 

w 

o 

< 

e 


Technical  Memorandum 


URUM-A  UNIFIED  RAD6ME 
LIMITATIONS  COMPUTER  PROGRAM, 

VOLUME  1- THEORETICAL 
BACKGROUND 


R.  K.  FRAZER 


qiF 


D D C 

fcEBgJME 

jj<  SEP  17  1978 


itSEUFE 

D 


THE  JOHNS  HOPKINS  UNIVERSITY  • APPLIED  PHYSICS  LABORATORY 


Approved  lor  public  distribution  unlimited. 


* *■-*  - - ' 


( 


APL/JHU 
TG  1293A 
JULY  1976 


I 

i 

i 


I ACCESSION  tor 

J 

NTIS 

Wi.lte 

OCli*  (3  1 

DOC 

lull  Sectlu  □ 1 

UKAMO'JHCEO 

JL'STIf  (CATION 

□ 

IT 

DISTRIOUTIOK/ATAILACIL 

ITT  CODES 

out, 

AVAIL.  *!»«/ 

,r  SPir.iri 

^ 'f 

I 

I 

I 

I 

I 

I 

I 

I 

I 


Technical  Memorandum 

DRUM -A  UNIFIED  RADOME 
LIMITATIONS  COMPUTER  PROGRAM, 
VOLUME  1- THEORETICAL 
BACKGROUND 


R.  K.  FRAZER 


D D C 

nr 


TPaEMSf 

i \ 

SEP  17  1976 

ilisissnrE 

jyj 

D - 


THE  JOHNS  HOPKINS  UNIVERSITY  ■ APPLIED  PHYSICS  LABORATORY 
Johns  Hopkins  Road,  Laurel,  Maryland  20810 


Approved  for  public  release,  distribution  unlimited. 


Unclassified 

SECURITY  CLASSIFICATION  OF  THIS  PAGE 

REPORT  DOCUMENTATION  PAGE 


2.  GOVT  ACCESSION  NO 


PLEASE  FOLD  BACK  IF  NOT  NEEDED 
FOR  BIBLIOGRAPHIC  PURPOSES 


3.  RECIPIENT'S  CATAIOG  NUMBER 


APL/ JHu/tG-1  29  3A-^(  ( - 


/ A 9."  I I M.  UP  IILIURI  . 

/URLIM  - A UNIFIED  RADOME  LIMITATIONS  COMPUTER  PROGRAM,  I ( j Technical/t^em:) 
JOLUME  lT^HEORETlCAL  JACKGROUND,  - * / V ' Ag  - i ) 

**  f 6 Koenouilft  /fa r 


6.  PFRFORMllT.  nRGlRFPORT  NUMBER 


t R.  K.Arazer 


m Nfl^7-72-C-44jBi 


9 PERFORMING  ORGANIZATION  NAME  & ADDRESS 

The  Johns  Hopkins  University  Applied  Physics  Laboratory 
Johns  Hopkins  Rd« 

Laurel,  MD  20810 

1 1 CONTROLLING  OFFICE  NAME  * ADDRESS 

Naval  Plant  Representative  Office 
Johns  Hopkins  Rd. 

Laurel,  MD  20810 

14.  MONITORING  AGENCY  NAME  a ADDRESS 

Naval  Plant  Representative  Office 
Johna  Hopkins  Rd. 

Laurel , MD  20810 


10.  PROGRAM  ELEMENT,  PROJECT.  TASK 
AREA  a WORK  UNIT  NUMBERS 


nu& 


IB.  SECURITY  CLASS.  ( of  thit  upon ) 


Unclassified 

IBs.  DECLASSIFICATION/DOWNGRADING 
SCHEDULE 


16  DISTRIBUTION  STATSMENT  [of  thii  fitport I 

Approved  for  public  release;  distribution  unlimited. 


17  OISTRIBU' ION  STATEMENT  [of  tho  obttroct  tnttrtd  in  Hock  X).  if  dlffortnt  from  fftpon) 


18  SUPPLEMENTARY  NOTES 


19.  KEY  WORDS  [Cononufon  rtvtrw  lido  if  noctuory  tnd  identify  by  btock  numbor I 

aerodynamic  heating  radoaes 

boresight  errors  thermal  stresses 

heat  transfer  trajectory  simulation 

missiles 

20  ABSTRACT  {Com.  nut  on  rtverss  sid*  if  noctuiry  tnd  identify  by  bfock  number) 

J^l'RLIM,  a unified  radome  limitations  computer  program,  has  been  developed  to  aid  the  radome  de- 
sign engineer  by  providing  a definition  of  the  maximum  flight  performance  capabilities  of  radome 
materials.  URLIM  numerically  determines  the  response  of  the  radome  to  aerodynamic  heating  and  load- 
ing. It  computes  the  following  as  functions  of  trajectory  time:  thermal  stress;  radar  boresight 

error  rates;  missile-radome  attachment  stresses  caused  by  maneuvers,  pressure,  and  drag  forces;  and 
the  onset  of  rodome  melting.  The  basic  output  of  the  program  is  a notation  of  trajectory  time  at 
which  the  radome  reaches  its  design  limitations.  Many  options  are  available  to  the  user  of  the 
URLIM  program  that  provide  a wide  variety  of  analysis  capability.  For  this  reason,  URLIM  may  also 
be  considered  as  a general  purpose  aerodynamic  heat  transfer  program  as  well  ae  a specific  purpose 
radome  limitations  program.  Volume  1 of  this  report  presents  the  theoretical  background  of  the 
analysis  techniques  used  in  URLIM;  Volume  2 provides  a detailed  explanation  of  how  to  use  URLIM. 


DD  1473 


63 / 6^0 


Unclassified 


SECURITY  CLASSIFICATION  OF  THIS  PA3E 


■■  j ZK-tsxa  


THE  JOHNS  HOPKINS  UNIVERSITY 

APPLIED  PHYSICS  LABORATORY 

LAUREL.  MARYUNO 


ABSTRACT 


URLIM,  a unified  radome  limitations  computer  program,  has 
been  developed  to  aid  the  radome  design  engineer  by  providing  a 
definition  of  the  maximum  flight  performance  capabilities  of  ra- 
dome materials.  URLIM  numerically  determines  the  response  or  the 
radome  to  aerodynamic  heating  and  loading.  It  computes  the  follow- 
ing as  functions  of  trajectory  time:  thermal  stress;  radar  bore- 

sight  error  rates;  missile-radome  attachment  stresses  causeO  by 
maneuvers,  pressure,  and  drag  forces;  and  the  onset  of  rado-.e 
melting.  The  basic  output  of  the  program  is  a notation  of  trajec- 
tory time  at  which  the  radome  reaches  its  design  limitatiors. 

Many  options  are  available  to  the  user  of  the  URLIM  progr-im  that 
provide  a wide  variety  of  analysis  capability.  For  this  it* ..son, 
URLIM  may  also  be  considered  as  a general  purpose  aerodynamic 
heat  transfer  program  as  well  as  a specific  purpose  radome  limita- 
tions program.  Volume  1 of  this  report  presents  the  theoretical 
background  of  the  analysis  techniques  used  in  URLIM;  Vo'  • me  2 
provides  a detailed  explanation  of  how  to  use  URLIM. 
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PREFACE 


The  purpose  of  this  volume  is  to  present  the  theoretical 
background  of  the  analysis  techniques  used  in  the  URLIM  program. 
This  discussion  will  consist  of  three  basic  parts:  first,  it  will 

present  the  mathematical  relationships  between  variables,  along 
with  their  assumptions  and  limitations;  second,  it  will  describe 
the  numerical  techniques  employed  by  the  program  to  solve  the  vari- 
ous equations;  third,  it  will  outline  the  data  management  tech- 
niques used  to  finW  the  various  solutions.  The  first  section  de- 
scribes the  physical  models  that  are  the  analytical  basis  of  the 
program.  As  such,  this  section  provides  a view  of  URLIM  that  will 
enable  the  prospective  user  to  assess  the  accuracy  and  completeness 
of  the  facilities  provided  in  the  program.  The  last  two  sections 
will  provide  an  overview  of  the  URLIM  program  organization  that  may 
give  the  interested  user  an  insight  to  the  generality  of  the  code. 
Moreover,  as  errors  will  inevitably  occur,  this  level  of  under- 
standing of  the  program's  organization  will  serve  as  an  aid  to  de- 
bugging. 
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1.  THEORETICAL  BASIS 


HEAT  TRANSFER  — CONDUCTION 

The  unified  radome  limitations  computer  program  (URLIM)  was 
developed  to  predict  heat  transfer  within  a solid.  Fundamentally, 
the  determination  of  the  time-dependent  temperature  field  within 
a closed  volume  will  depend  on  the  rate  of  heat  conducted  orthog- 
onal’./ across  the  surface  of  the  volume  and  the  rate  of  change  of 
internal  energy  within  the  volume.  V-’e  will  assume  the  arbitrary 
case  of  a volume  V with  some  temperature  existing  at  every  point 
within  V.  In  the  absence  of  any  effects  from  electric  or  magnetic 
fields,  surface  tension,  or  chemical  reactions,  the  energy  in  a 
solid  is  a function  of  its  temperature  such  that 


du 

dt 


3 (pc  T) 

E 

3t 


(1) 


where 

du/dt  is  the  time  rate  of  change  of  internal  energy  in  a 
differential  volume,  dV, 

p is  the  density  of  the  material, 

Cp  is  the  specific  heat  of  the  material, 

T is  the  magnitude  of  the  temperature  field  at  the  dif- 
ferential volume  dV  (i.e.,  the  temperature  of  dV) , 
and 

t is  the  independent  variable  time. 

The  movement  of  heat  through  any  point  in  the  volume  is  related  to 
the  gradient  of  the  temperature  field  within  the  volume;  that  is, 


Q - -k7T 


(2) 


where 


Q is  the  heat  flow  vector. 
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7T  is  the  temperature  field  gradient  vector,  and 

k is  the  proportionality  factor  that  can  only  be  a func- 
tion of  temperature  (i.e. , k ■ f(T)). 

Equation  (2)  states  that  heat  flows  within  a solid  in  the  opposite 
direction  to  the  temperature  gradient  and  that  the  magnitude  of 
that  flow  is  related  by  the  material  property  k,  which  may  be 
temperature  dependent . 

With  Equation  (2)  describing  the  flowjof  heat  at  any  point 
in  the  volume,  the  divergence  of  the  vector  Q is  written  as 


7 • Q ■ V * k7T 


(3) 


and  is  thought  of  is  the  change  in  magnitude  cf  heat  flow  with  re- 
spect to  any  point  in  the  volume.  In  a conservative  system  (as 
has  been  assumed  here)  this  spatial  change  in  heat  flow  at  a point 
in  the  volume  is  equal  to  the  temporal  change  in  internal  energy 
of  that  point.  Using  2qs.  (1)  and  (3)  the  following  equation  rep- 
resents this  8tate«ucr.L  of  the  conservation  of  energy: 


7 • kTT 


3(pcpT) 

3t 


(4a) 


The  material  property  product  (pcp)  is  allowed  to  depend  explicitly 
only  on  temperature  (i.e.,  pcp  » f(T)),  so  that 

7 • k7T  - -pc  . (4b) 

p at 


The  operator  7 is  the  vector  gradient  operator  and  must  be 
represented  in  an  orthogonal  coordinate  system.  It  can  have  the 
following  definitions. 

Cartesian  Coordinates: 


7 


(5) 
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where  x,  y,  and  z are  the  unit  vectors  in  each  of  the  coordinate 
directions. 

Cylindrical  Coordinates: 


V 


9 


where 

r ia  the  radial  coordinate, 

6 is  the  circumferential  coordinate,  and 
z is  the  axial  coordinate. 

Spherical  Coordinates: 


(6) 


V - 


3_ 

Sr 


+ !»-*  + —I—  i_ 

r 3<  ’ r sin  ♦ 30 


e 


(7) 


where 

r is  the  radial  coordinate, 

$ is  the  co-latitudinal  coordinate,  and 
8 is  the  longitudinal  coordinate. 

Equation  (4)  is  the  governing  relationship  between  the 
spatial  ’»nd  temporal  changes  of  temperature  within  the  body.  Spe- 
cific soj.jt.ion8  to  the  sets  of  differential  equations  Indicated  in 
Eq,  (4)  depend  on  the  specification  of  the  conditions  that  exist 
at  the  boundary  surface  of  the  volume  and  whether  or  not  there  are 
independent  heat  sources  or  sinks  present  within  the  volume.  Sec- 
tion 2 of  this  volume  will  show  how  the  general  relationships  ex- 
pressed here  are  modeled  and  the  temperature  field  solved  for  as  a 
function  of  time. 


HEAT  TRANSFER  — CONVECTION 

In  situations  where  the  volume  being  considered  is  subjected 
to  a fluid  flowing  at  its  surface,  heat  transfer  will  occur  across 
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the  surface  due  to  conduction  within  the  fluid,  viscous  friction 
dissipation  within  the  fluid  near  the  surface,  and  molecular  mo- 
tion in  the  boundary  layer  (i.e.,  convection).  In  general,  the 
static  temperature  of  the  fluid  and  the  surface  temperature  of  the 
solid  will  be  different  and  there  will  be  a potential  for  moving 
heat  across  the  surface.  Figure  1 shows  a schematic  diagram  of  a 
fluid  flowing  past  a solid  surface  with  an  appropriate  coordinate 
system  and  with  the  important  parameters  listed.  If  the  bulk  of 
the  fluid  stream  in  the  vicinity  of  the  surface  is  hotter  than  the 
surface,  there  will  be  a spatial  variance  of  temperature  through 
the  stream  as  showr  in  the  figure.  Due  to  the  viscous  effects  of 
the  fluid,  the  fluid  particles  at  the  wall  surface  will  be  assumed 
stationary  with  respect  to  the  wall.  The  velocity  and  temperature 
variances  as  one  travels  away  from  the  well  are  as  shown  in  the 
figure.  The  points  at  which  the  velocity  and  temperature  of  the 
gas  reach  99 X of  their  freestream  values  describe  the  boundary 
layer  edge;  the  two  points  (one  for  temperature  and  one  for  veloc- 
ity) will  be  asst  '•d  to  coincide  spatially  at  every  point  along 
the  body.  The  maximum  boundary  layer  temperature  is  called  the 
"recovery"  or  "adiabatic  wall"  temperature  and  is  defined  as  the 
fluid  temperature  where  the  derivative  with  respect  to  the  local 
wull  normal  (3Tf/3y)  is  zero. 

Predictions  of  the  aerodynamic  heating  as  described  so  far 
involve  calculating  the  recovery  temperature  and  estimating  the 
proper  dependence  of  the  convective  flux  (q)  on  the  recovery  tem- 
perature. In  the  discussions  that  follow,  the  basic  relation 
shown  below  will  be  used  to  define  the  heat  flux  per  unit  area  to 
the  solid  surface: 

% m hl(1r  ' V • <8) 


where 

h^  is  a heat  transfer  coefficient  that  will  relate  heat 
1 flow  to  a difference  in  enthalpy, 

i is  the  enthalpy  of  the  fluid  at  the  recovery  tempera- 
r ture,  and 

i is  the  enthalpy  of  the  fluid  at  the  wall  temperature. 
Enthalpy,  a valid  material  property,  is  defined  as 
i * u + pv, 

where  u is  the  internal  energy,  p the  pressure,  and  v the 
specific  volume  (reciprocal  of  density).  Enthalpy  is,  in  general, 
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Freestream  condition 


Density  - 
Viscosity  - Moo 
Temperature  - T^, 
Pressure  - 
Velocity  - Vqq 


Shock  wave  for  compressible  supersonic  flows; 
implies  local  fluid  properties  that  will  depend 
^ on  coordinates  X,  V and  body  shape 

(See  detail  below,  left) 

External  surface  of  body 

First  internal  nodes 


Body  coordinates 


Local  flow  condition  * Velocity  (Vj) 

» Temperature  (T|) 


Local  coordinates 


Boundary  layer  edge 
Heated  surface  at  Tw 


Boundary 


Velocity 


Two  thermal  nodes  with  similar  local  flow 
at  boundary  layer  edge 


Boundary 
layer  edge 


For  a heated  wall 


For  an  adiabatic  wall 


Tw  Tj  ’recovery 


Fluid  temperature,  Tf 


Fig.  1 Terminology  for  Aerodynamic  Heating  Boundary  Conditions 
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dependent  on  temperature  and  pressure.  We  will  assume  knowledge 
of  the  fluid  properties,  Including  enthalpy  versus  fluid  tempera- 
ture (Tf ) and  fluid  pressure  (Pf) , so  that  it  remains  to  calcu- 
late hi  and  the  recovery  enthalpy.  In  the  following  discussion 
the  fluid  will  be  assumed  to  be  a gas  (usually  air) . 

With  a compressible  fluid  as  the  medium  there  is  the  possi- 
bility that  the  freestream  velocity  is  larger  than  the  acoustic 
velocity  in  the  gas  and  a shock  wave  may  be  present  in  the  vicinity 
of  the  wall  surface  being  considered.  If  this  is  the  case,  the 
local  boundary  layer  edge  velocity  (V^)  must  first  be  calculated. 
This  calculation  involves  a priori  taowledge  of  the  ratio  of  pres- 
sures across  the  shock  wave  as  well  as  the  ratio  of  Mach  numbers. 
(The  Mach  number  (M)  is  the  ratio  of  actual  gas  velocity  to  the 
acoustic  velocity  of  the  gas  at  any  point.)  Further,  the  distance 
from  the  shock  wave  to  the  local  body  position  and  the  shape  of 
the  body  have  effects  on  the  Mach  and  pressure  ratios  mentioned 


In  other  v-  Is, 

Vp0  * f (Mq,x)  and 

(9) 

Ml/M0  “ 8(Mq,x)  , 

(10) 

with  x indicating  the  position  along  the  body,  1 indicating  the 
local  condition,  0 indicating  the  freestream  condition,  and  the 
other  items  as  shown  in  Fig.  1.  Since  energy  is  conserved  across 
a shock  wave,  the  following  equation  will  hold: 


i 


1 


(11) 


which  simply  sun>e  the  sensible  and  kinetic  energies  on  either  side 
of  the  ahcck  wave.  If  the  gas  is  assumed  to  be  ideal  except  for  a 
compressibility  factor  (Z)  defined  as 

Z « pv/RT, 

then  the  sonic  velocity  (a)  can  be  shown  to  be 

a - y}y$zkT  , (12) 

where 

Y is  the  ratio  of  specific  heats  (cp/cv)  for  the  gas, 
which  are  defined  as 
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c ■ di/dT  and 
P 

c * du/dT  • 
v 

The  gas  properties  y and  Z will  be  assumed  to  be  known  as  func- 
tions of  temperature  and  pressure.  Using  the  above  oof initions 
Eq.  (11)  is  written  as 


*1  + 


VlRTlMl 

2 


10  + 


X 

2 


(13) 


With  knowledge  of  the  freestream  conditions  of  temperature  and 
pressure,  the  gas  properties  i,  y,  Z,  and  the  functions  indicated 
by  Eqs.  (9)  and  (10),  Eq.  (13)  is  seen  as  an  implicit  equation  in 

Tl- 


When  Eq,  (13)  is  solved  for  T].,  the  local  flow  conditions 
are  then  fully  defined  and  the  recovery  condition  is  next  con- 
sidered. The  recovery  temperature  within  the  boundary  layer  must 
lie  somewhere  between  the  total  temperature  and  the  local  boundary 
layer  edge  temperature.  Thinking  in  terms  of  enthalpy  instead  of 
temperature  this  notion  is  expressed  as 


i 


r 


(14) 


where  r is  the  recovery  factor,  defined  as 


r 


(15a) 


To  determine  r it  is  necessary  to  resort  to  empirical  correla- 
tions of  measured  data.  Reference  1 states  the  following  relation- 
ship for  r: 


r « \ Pr*  for  laminar  flow  (15b) 


Ref.  1.  R.  E.  Wilson,  Handbook  of  Supersonic  Aerodynamics, 
Sections  13  and  14,  "Viscosity"  and  "Heat  Transfer  Effects," 
NAVORD  Report  1488,  Naval  Ordnance  Laboratory,  White  Oak,  MD, 
August  1966. 
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r * yj  f r*  for  turbulent  flow,  (15c) 

where  ?r*  is  the  Prandtl  number  of  the  airstream  evaluated  at  the 
reference  temperature  (reference  denoted  by  *)  and  the  Pr  number  is 


Pr  - cp  u/k  , 

p being  the  viscosity.  We  have  now  traded  one  unknown  for  another; 
that  is,  the  recovery  factor  is  defined  In  terms  of  another  un- 
known condition  — the  reference  temperature  (T*).  The  reference 
temperature  is  a fictitious  boundary  layer  temperature  that  is  used 
to  evaluate  the  physical  properties  of  the  stream.  Eckert  (Ref.  2) 
proposes  the  following  relationship  for  determining  the  reference 
condition  in  terms  of  enthalpies: 

i*  - + t1  + 0.22  r Vx2)/2  , (16) 


which  is  seen  to  be  an  average  of  enthalpies  and  a kinetic  energy 
factor.  Equation  (16)  now  involves  the  recovery  factor  so  an  ex- 
plicit solution  for  r from  Eqs.  (14),  (15),  and  (16)  cannot  be 
made;  the  relationships  must,  be  solved  simultaneously. 

The  results  uf  solving  Eqs.  (14)  through  (16)  are  the 
reference  enthalpy  (i*)  and  the  recovery  enthalpy  (ir) . By  use 
of  the  fluid  property  tables  the  reference  and  recovery  enthalpies 
can  be  used  to  obtain  the  reference  and  recovery  temperatures  (T* 
and  Tr,  respectively).  At  this  point  we  can  turn  to  heat  trans- 
fer correlations  to  evaluate  the  h^  of  Eq.  (8).  The  Stanton  num- 
ber is  defined  as 


St  - 


»*vi(1r  - V 


(17) 


Ref.  2,  E.  R.  G.  Eckert,  "Survey  of  Boundary  Layer  Heat 
Transfer  at  High  Velocities  and  High  Temperatures,"  Technical  Re- 
port 59-624,  Wright  Air  Development  Center,  Dayton,  OH,  April  1960. 
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and  Colburn  (Ref.  3)  notes  that 

t 

St  - cf  (Pr*)0  . 


(18) 


The  term  cf  is  the  friction  factor  and  is  further  defined  as 


w 


pVl  /28c 


(19) 


where 

t is  the  stream  shear  stress  at  the  wall  surface, 
w 

gc  is  32,174  ft-lbm/lbf-s2,  and 

0*  is  an  empirical  constant. 

Experimental  correlations  (Refs.  2 and  3)  have  shown  that 

cf  - C Rea  , (20) 


where  C and  a are  again  empirical  constants,  and  Re  is  the 
Reynolds  number  (pVjx/p).  Combining  Eqs.  (17),  (18),  and  (20) 
gives 


P*VT 

44 


/ 4 

C Rea  Pr*& 


(21a) 


Using  the  definitions  for  Re  and  Pr,  and  the  following  for  Nusselt 
numbers 


Nu  * (c  x h )/k  , 
P i 


Ref.  3.  A.  P.  Colburn,  "A  Method  of  Correlating  Forced  Con' 
vection  Heat  Transfer  Data  and  a Comparison  with  Fluid  Friction," 
Trans.  AIChE,  Vol.  29,  1933,  pp.  174-210. 
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the  general  relationship 


Nu  « C Rea  PrB 


(21b) 


can  be  written  where  the  constants  a and  8 are  different  from 
a and  8*  and  are  given  values  according  to  Table  1. 

Given  the  position  on  the  body  (x) ; the  previously  calcu- 
lated values  for  Tr,  7*,  and  Vi;  and  the  stream  property  func- 
tions; Eq.  (21)  or  (21a)  is  solved  for  h^  and  Eq.  (8)  is  used  to 
find  the  heat  flux  to  the  wall  (qw) . 

In  summary,  the  following  steps  are  required  to  calculate 

V 


1.  Know  i,  Z,  y>  Pr,  and  p for  the  stream  fluid  versus 
T and  P ; 

2.  Know  the  Mach  and  pressure  ratios  (Eqs.  (9)  and  (10)) 
between  the  freestream  and  the  local  boundary  layer 
edge; 

3.  Calculate  the  local  boundary  layer  edge  static  tempera- 
ture (Tj)  from  Eq.  (13); 

4.  Solve  for  the  recovery  factor  (r),  recovery  temperature 
(Tr) , and  reference  condition  (T*)  (Eqs.  (14),  (15), 
(16)); 

5.  Using  properties  at  the  reference  temperature,  solve  for 
hi  using  Eq.  (21);  and 

6.  Substitute  hj,  ir,  and  iw  in  Eq.  (8)  to  solve  for 

Whether  the  laminar  or  turbulent  coefficients  in  Table  1 are  used 
depends  generally  on  the  value  of  the  Reynolds  number;  low  Re  im- 
plies laminar  flow  and  high  Re  implies  turbulent.  The  transition 
value  of  Re  must  be  decided  a priori. 

In  the  case  of  heating  of  a surface  that  is  normal  to  the 
freestream  (i.e.,  stagnation  points)  the  methods  of  Squire  and 
Sibulkin  are  used  (Refs,  6 and  7).  In  these  techniques  it  is 

Ref.  6.  S.  Goldstein,  Modern  Developments  in  Fluid  Dynamics, 
First  Edition,  Vol.  2,  Oxford  Univ.  Press,  London,  1938,  p.  631. 

Ref.  7.  M.  Sibulkin,  "Heat  Transfer  near  the  Forward  Stag- 
nation Point  of  a Blunt  Body,"  J.  Aeronautical  Sciences,  August 
1.952.  ” 
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possible  to  show  that  if  ti, : square  root  of  the  stream  velocity 
gradient  at  the  stagnation  point  is  'incorporated  in  the  constant 
C (Table  1)  then  proper  correlation  with  measured  data  is  ob- 
tained. In  these  stagnation  point  heating  cases,  it  was  found 
that  the  term 


v'd/v. 


was  required  in  the  constant  C of  Eq.  (21)  where 

V*  is  (3V/3y)|y  » 0,  or  the  velocity  gradient  of  the 
stream  evaluated  at  the  wall  (c.f.  Fig.  1), 

D is  the  characteristic  diameter  of  the  body,  and 

V is  th<  freestream  velocity. 

00 

The  dependence  of  v'  on  Mach  number  has  been  derived  for  a real 
gas  using  Newtonian  flow  assumptions.  This  plot  is  given  in 
Fig.  2 and  will  be  required  for  use  when  stagnation  point  heating 
is  being  considered. 

The  aerodynamic  heat  transfer  relationships  discussed  above 
are  solved  in  the  SAERO  routine  where  the  final  result  is  simply 
the  temperature  for  the  surface  node  in  question.  The  temperature 
is  solved  for  by  making  an  energy  balance  at  the  heated  surface 
per  unit  area  as  follows  (c.f.  Fig.  1): 


« oe 

w 


+ 4 


con 


t 


where ; 

q , is  an  independently  specified  heat  flux  to  the 
ra  surface, 

T^  is  the  wall  temperature, 

q is  the  heat  conducted  from  the  surface  into  the 

con  material  below  the  surface, 

o is  the  Stephan-Boltzmann  constant  (see  the  follow- 

ing section  on  radiation  heat  transfer),  and 

ew  is  the  total  normal  emissivity  of  the  wall  surface. 
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In  this  relationship,  qw  (the  aerodynamic  heat  transfer  to  the 
wall)  is  explicitly  related  to  Tw  through  the  enthalpy  iw  (c.f. 

Eq.  (8)).  qra(j  Is  independently  specified  and  is  therefore  a con- 
stant at  any  instant  in  time.  Tw^  is  an  obvious  function  of  Tw» 
and  qcon  depends  on  Tw  as  follows: 


l ■ k (T  - T. )/L  . 

con  w w i w-i 


where  T*  is  the  inner  wall  temperature  within  a node  at  a distance 
Ly.j  from  the  surface.  The  energy  balance  equation  above  is  then 
seen  as  an  implicit  equation  in  Tw,  if  it  is  assumed  that  Ti  is 
known.  At  each  time  increment  in  the  transient  solution  technique 
the  temperature  T-t  from  the  previous  time  point  is  used  so  that  the 
conduction  equation  just  above  can  be  solved  for  qCon  anc*  the  en- 
ergy balance  equation  can  be  solved  for  Tw. 

At  thi'.  point  we  will  also  note  that  the  solution  for  r, 

Tr.  and  T*  (Eqs.  (14),  (15«),  and  (16))  is  made  an  explicit  rela- 
tion in  the  same  way  as  the  internal  conduction  equation.  That 
is,  the  term  r (Eq.  (15a))  is  used  from  a previous  time  step  so 
that  ir  (Eq.  (14))  and  i*  (Eq.  (16))  are  found  without  the  need 
for  a simultaneous  solution  of  the  three  equations. 


HEAT  TRANSFER  - RADIATION 

Every  substance  at  any  temperature  above  absolute  zero 
emits  electromagnetic  (EM)  wave  energy  that  will  cause  heat  to  be 
absorbed  by  any  other  surface  upon  which  the  energy  impinges. 

Heat  is,  on  an  atomic  level,  a measure  of  the  intensity  of  molecu- 
lar vibrations.  The  molecules  in  a solid  are  groups  of  positive 
and  negative  electrical  charges  so  their  vibration  at  a surface 
will  cause  EM  radiation.  Conversely,  EM  waves  of  the  proper  fre- 
quency impinging  on  a surface  will  excite  the  molecular  vibrations 
of  that  surface.  In  1879  Stephan  observed  the  following  relation- 
ship for  the  energy  flux  radiated  by  a body  at  temperature  T with 
a perfect  surface  (that  is,  a surface  that  emits  the  EM  radiation 
without  loss) : 


Eb‘ 


oT 


(22a) 


Equation  (22a)  was  also  theoretically  derived  by  Boltzmann  and  is 
termed  the  Stephan-Boltzmann  law;  o is  the  Stephan-Boltzmann 
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constant.  For  real  surfaces,  the  amount  of  energy  radiated  is 
less  than  the  amount  given  by  Eq.  (22a)  and  an  emissivity  factor 
(e)  is  generally  incorporated: 

E - eoT4  . (22b) 

g 


The  subscripts  b and  g in  Eqs.  (22a)  and  (22b),  respectively, 
represent  "black-body"  (perfect)  and  "grey-body"  (imperfect)  sur- 
faces, respectively.  Because  solid  bodies  are  aggregates  of 
molecules  that  are  subject  to  the  laws  of  quantum  physics  and 
quantum  probabilities,  the  energy  emitted  at  a particular  tempera- 
ture will  be  distributed  over  a wide  range  of  wavelengths;  l.e., 
the  surface  molecules  vibrate  over  a range  of  frequencies.  In 
1900  Plank  derived  the  following  relation  for  the  energy  emitted 
at  a specific  wavelength  by  a body  at  a temperature  T (c.f. 

Chapter  5 of  Ref.  ?.) 


V 


-5 


exp(C2/AT)  - 1 


(23) 


where  and  C2  are  constants.  Equation  (22a)  is  then  seen  as  an 
integration  of  Eq.  (23)  over  all  wavelengths;  i.e., 


E 


b 


J EbXdX  " C1  f exp(C2/AT)  - 1 
0 0 


The  implication  of  Eq.  (23)  is  that  for  real  surfaces  the 
emissivity  may,  in  general,  exhibit  a dependence  on  A and  T 
since  a surface  could  be  more  emissive  at  one  wavelength  than  an- 
other. 


For  a surface  exposed  to  the  EM  radiation  from  another  hot 
body  the  fraction  absorbed  (a)  is  related  to  the  fraction  reflected 
(r)  and  transmitted  (x)  by  a simple  expression  of  the  conservation 
of  energy: 


Ref.  8.  F.  Kreith,  Principles  of  Heat  Transfer,  Second  Edi- 
tion, Section  5,  International  Textbook  Co.,  Scranton,  PA,  1965. 
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a + T + «.  * 1 


where 


a is  the  absorptivity, 

t is  the  transmissivity,  and 

r is  the  reflectivity, 

all  of  which  may  depend  on  T and  X. 

In  practical  situations  where  radiation  heat  transfer  is 
significant,  the  range  of  temperatures  are  such  that  the  wave- 
lengths of  the  EM  energy  involved  are  roughly  between  0.5  * 10”® 
and  20  x 10”®  m.  This  EM  waveband  is  termed  the  infrared  (IR) 
bard.  A simplifying  assumption  that  will  be  made  here  is  that  all 
real  materials  hn  values  of  r,  e,  and  a that  are  constant 
over  the  IR  band  and  can  depend  only  on  the  temperature  of  the 
surface.  Further,  the  so  called  grey-body  assumption  will  be 
made;  that  is,  at  a given  temperature  the  emissivity  of  a material 
will  be  equal  to  the  absorptivity.  Given  these  assumptions,  the 
net  exchange  of  heat  energy  between  the  two  surfaces  shown  below 
can  be  described. 


It  is  desired  to  know  the  net  heat  flux  per  unit  area 
^net^  across  the  imaginary  plane  spaced  between  the  two  given 
planes.  This  flux  will  be  the  total  radiated  heat  per  unit  area 
from  the  surface  at  T^  less  the  heat  radiated  from  the  surface  at 
T2  plus  a consideration  for  the  energy  reflected  at  each  surface. 
In  algebraic  terms: 
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q ■ eoT  ^ - eoT,.4  + (1  - «)eoT  4 - 
net  1 l i 

(1  - a)eaT  4 + (1  - a)2  eoT^  - (1  - cx)2ta1^  + ..  . (24) 


The  first  two  terms  are  the  primary  radiation  terms,  the  third 
term  represents  the  heat  radiated  by  surface  1 but  reflected  at 
surface  1,  the  fourth  term  is  the  heat  similarly  reflected  from 
surface  2,  and  the  successive  terms  represent  the  continuing  re- 
flections that  proceed  indefinitely.  Equation  (24)  can  be  simpli- 
fied to: 


’net  ‘ “ 


Tj4  - T24  + Tj4  ^ <1  - «)“  (-l)” 

n*l 


- t24  2 U ^ ("1>n 


n*l 


(25a) 


or 


q ■ eo 
^net 


(Tl4  - T24)S  (1  ' ('1>n 

n"0 


(25b) 


The  infinite  sum  is  equal  to: 


£ a - a>n  <-i,n  * 2^  • 


n-0 


(26) 


Using  the  assumption  that  a • e,  and 


Tl4  " T24  ^ 4 t3  (T1  ‘ V * 
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then  Eq.  (25a)  becomes 


i o — 521 — (x  - x ) 
’net  (2  - e)  U1  2}  ’ 


(27a) 


where  T3  - (^  + T2)  ^ + T £ J 


If  we  now  consider  the  distance  between  the  surfaces  (d) 
and  rewrite  Eq.  (27a)  as 


’net 


eoT3d 
(2  - e) 


(T1  - T2) 


(27b) 


then  the  term  in  the  square  brackets  can  be  regarded  as  the  "ef- 
fective" thermal  conductivity  of  the  space  (d)  between  1 and  2; 
i.e. , 


AT 


*net  eff  d 


with 


eoT3d 
Ceff  " (2  - e) 


(27c) 


The  above  treatment  is  for  parallel,  infinite  plates  but 
similar  analyses  can  be  done  considering  the  geometrical  differ- 
ences for  arbitrary  surfaces  and  the  result  is  to  add  another  fac- 
tor called  the  view  factor  (F)  to  Eq.  (27b);  ie., 


eoT3F 
eff  " (2  - e) 


(27d) 


Values  of  F have  been  tabulated  for  various  conventional  surface 
arrangements;  Ref.  8 gives  some  typical  values. 
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While  Eq.  (27d)  gives  an  effective  conductivity  for  a 
linearized  temperature  function,  Eq.  (25b)  could  just  as  well  be 
written 


q 


net 


» 


with 


eeff  " £/(2  " e)‘ 


In  Section  2 of  this  volume  the  solutions  to  these  rela- 
tionships will  depend  on  having  the  values  eeff  mentioned  above 
supplied  a priori. 


THERMAL  STRESSES 

Subroutines  in  URLIM  can  calculate  thermal  stresses  in 
arbitrarily  defined  cylindrical  wall  sections.  The  method  em- 
ployed is  that  of  Rivello  (Ref.  9)  and  is  briefly  described  here. 

A cylinder  of  infinite  length  subdivided  into  n concentric 
cylindrical  sections  i3  shown  schematically  in  Fig.  3.  Also  shown 
in  Fig.  3 Is  a definition  of  the  terms  used  in  this  discussion. 

The  following  assumptions  are  made: 

1.  The  elastic  modulus  relating  the  stress  in  the  material 
Lo  the  strain  is  constant  within  each  region. 

2.  Poisson's  ratio  is  constant  throughout  all  regions. 

3.  The  materials  that  may  make  up  the  cylinder  are  iso- 
tropic. 

4.  The  cylinder  is  restrained  from  motion  in  the  axial 
direction;  i.e.,  the  axial  strain  is  zero  at  all 
values  of  r. 

5.  The  radial  distribution  of  temperature  is  known  at  all 
values  of  radius. 


Ref.  9.  R.  M.  Rivello,  "Thermal  Stress  Analysis  of  Sand- 
wich Cylinders,"  APL/JHU  TG  721,  August  1965. 
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NOTES; 

•Total  wall  thickness  (bn+^  — b-j ) is  divided  into  n arbitrarily  sized  regions, 
creating  n+1  interfaces. 

•Temperatures  are  provided  at  m arbitrarily  spaced  radius  values  and  arc-  constant 
with  respect  to  the  axial  and  circumferential  directions. 

• Radial  displacements  from  an  initial  isothermal  state  occur  for  ea-h  interface,  uj  • 

• Mechanical  propert  ies  are  assumed  constant  within  regions  at  values  that  are 
temperature  averaged  with  respect  to  the  radial  temperature  variance. 


Fig.  3 Thermal  Stress  Geometry  and  Nomenclature 
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If  T(r)  is  the  temperature  at  radius  r and  Tq  is  the  uniform  tem- 
perature for  the  cylinder  at  which  no  stresses  are  present  then 
T ■ T(r)  - Tq  will  define  the  variable  T in  the  following  devel- 
opment . 


The  radial  displacement  of  any  point  within  the  ith  region 
of  the  cylinder  in  Fig.  3 due  to  the  temperature  gradient  T is 
(from  Ref.  10) 


u 


i 


1 + v 1 

1 - v r 


(aT)r  dr  + Ci  L 


(28) 


The  radial  (r),  tangential  (6),  and  axial  (z)  stresses  are  given 
(respectively)  as 


(aT)r  dr  + 


1 + v 


(29) 


\ i f (otT)Ei 

ofl  - / (aT)r  dr  - - 

6.  1 - v 2 / v ' 1 - v 

i r .J 


1 + vl  1 - 2v  r2  J 


, and 


(«T)E1  2vEiCi,i 
°zi  1 ” v + (1  + v)(l  - 2v) 


(30) 

(31) 


In  these  equations  Cj  ^ and  2 are  constants  of  integration  that 
apply  to  the  ith  region. 


Ref,  10.  S,  Timoshenko  and  J.  N.  Goodier,  Theory  of  Elastic- 
ity. McGraw-Hill,  New  York,  i951. 
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The  coefficient  of  expansion  has  been  brought  under  the 
integral  sign  in  writing  these  equations  since  it  is  always  per- 
missible to  consider  oT  (the  unrestrained  strain  due  to  tempera- 
ture) as  a single  quantity.  As  a result  a need  not  be  constant 
over  the  range  of  temperature  changes. 

Since  there  are  n regions,  a total  of  2n  constants  of  inte 
gration  must  be  found.  These  constants  are  determined  so  that 
equilibrium  and  compatibility  are  satisfied  on  the  faces  and  inter 
faces  of  the  layers.  Equilibrium  on  the  inner  and  outer  faces 
requires  that 


a 


r 


and 


a 

r 


Pn+1 


(32) 


(33) 


where  pj_  and  pn+1  are  the  pressures  on  the  inner  and  outer  faces. 
At  each  of  the  interfaces  the  equilibrium  condition, 


i » 1 to  (n  - 1)  , 


(34) 


and  the  compatibility  condition , 


Ullb1+1  ’ Vllw  ' 


i ■ 1 to  (n  - 1)  , 


(35) 


must  be  satisfied,  Substituting  Eqs.  (28)  and  (29)  into  Eqs.  (32) 
through  (35)  gives 


P,(l  + v) 


( . i—  \ c -JLC  - 

\1  - 2v  / 1,1  b 2 ^1,2  E1 


(36) 
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i 


i 

l 

l 

i 

i 

I 

1 

l 


I 

1 

1 

I 

I 

I 

S 


r-S  Cn.2  ’ » + *> 

‘W 


_ fn+1 

(b., \ 2 En 

TiTxy 


Ji+1 


c,.,  , + 


,h  \2  1,2  E.  (1  - 2v)  i+1,1 

(bi+l)  1 


(37) 


Ei+1 

■ At  — o C. , . , - (1  + v) --- -r  ; i-1,2 (n-l),  (38) 

il(b1+1>2  1+1 ■ 2 <b1+l>2 


and 


1+1.1  (h  ^ 

^i+l* 


- - (1  + v) 


i*l 1 2 1 • • • t (n-l) 


where 


(aT)r  dr  . 


(39) 


(AO) 


Solution  of  the  2r.  simultaneous  equations  generated  by 
Bqs.  (36)  through  (40)  gives  the  ^ *nc*  Cj  o coefficients.  Sub- 
stitution of  these  into  Eqs.  (29)  through  (*i)  yields  or.,  o^  and 
oz^  for  restrained  ends. 

The  resultant  axial  force  for  restrained  ends  is 
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which,  by  using  Eq.  (31),  becomes 


R - — 2ir 


i & /1+1 


2vEiCi  1 
(aT)r  dr i-ixk. 


(1  + v)(l  - 2v) 


r 


r dr 


From  Eq.  (AO)  we  find 


vE.C, 


r . _2  V E a - 
it  2~i  i En  a + v)d  - 


2v) 


<bi+l)2  * bi2 


i-1 


(Al) 


To  determine  the  axial  stresses  for  unrestrained  ends  we  determine 
the  stresses  due  to  an  axial  force  of  -R  and  superimpose  these  upon 
the  previously  determined  stresses.  The  stresses  due  to  -R  can  be 
shown  to  be 


In  the  ith  region  the  radial  and  circumferential  stresses  for  un- 
restrained ends  are  then  given  by  Eqs.  (29)  and  (30).  The  axial 
stress  is 


ctTE 

rr?  + 


2vEici.i 


(1  + v)(l  - 2v) 


Ei(R/ir) 

<Vi)2-V*] 


. (A2) 


The  thermal  stresses  as  developed  here  are  used  in  the 
URLIM  program  to  evaluate  the  thermal  stress  failure  levels  of 
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radomes.  An  obvious  further  assumption  is  that  the  region  of 
critical  stress  in  the  radome  behaves  like  the  cylinder  modeled 
here.  Studies  by  others  (Ref.  11)  have  shown  that  if  the  point 
of  interest  on  the  radome  is  sufficiently  far  from  the  nose  re- 
gion, then  the  cylindrical  assumption  is  sufficiently  accurate 
(i.e.,  less  than  10%  error).  It  is  worth  noting  that  the  radial 
dimensions  that  compare  between  the  cylinder  and  the  radome  should 
be  made  along  the  local  normal  to  the  radome  profile  and  not  along 
the  radius  of  revolution  of  the  radome  profile. 


BORESIGHT  ERRORS 

The  prediction  of  boresight  error  rates  for  radar  transmis- 
sion through  streamlined  radomes  is  regrettably  imprecise.  This 
situation  is  due  in  large  part  to  the  number  of  significant  varia- 
bles involved,  such  as  radar  frequency,  antenna  design,  antenna 
placement  within  the  radome  and  incidence  angle  variations,  ra- 
dome shape,  the  magnitude  and  temperature  dependence  of  the  radome 
dielectri"  constant,  and  the  transmission  loss  of  the  material. 
Theoretical  prediction  techniques  that  are  sufficiently  general 
and  complete  are  consequently  quite  complicated.  To  add  to  the 
confusion,  the  prediction  techniques  often  do  not  correlate  well 
with  observed  data  (c.f.  Section  8 of  Ref.  12).  Therefore,  the 
philosophy  of  the  radome  error  prediction  method  used  here  is 
based  on  experimental  radome  data.  Interpolation  of  the  avail- 
able data  is  made  to  determine  the  empirical  equation  that  de- 
scribes the  effects  of  antenna  aperture  size,  dielectric  constant, 
wall  thicknese,  and  wavelength  changes  on  a theoretically  perfect 
room-temperature  radome  design.  An  attempt  is  also  made  to  in- 
clude the  effects  of  missile  dynamics  on  error  slope  requirements 
by  averaging  the  experimental  radome  error  over  a fixed  20°  gim- 
bal  period.  Before  a discussion  of  the  boresight  error  analysis 
is  given,  a brief  account  of  the  radar  power  loss  due  to  trans- 
mission through  a radome  will  be  given. 

Transmission  Losses 


A radome  must  pass  electromagnetic  radiation  efficiently. 

A reduction  of  radome  transmission  efficiency  degrades  the  overall 
range  product  or  effective  noise  figure  of  an  airborne  radar. 
Definition  of  a radome  limit  from  transmission  loss  requires 


Ref.  11.  R.  M.  Rivello,  Comparisons  of  Radome  Stress  Solu- 
tions," APL/JHU  EM-3989,  August  1965. 

Ref.  12.  L.  B.  Weckesser  et  al.,  "Environmental  Limitations 
of  Alumina,  Fused  Silica,  and  Pyroceram  9606  Radomes,"  APL/JHU 
TG  865,  May  1967. 
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knowledge  of  target  size  and  miasile-target  space  geometries. 

The  limitation  of  a radome  from  a transmission  loss  analysis  is 
very  vague  since  transmission  losses  simply  reduce  acquisition 
range.  It  is  conceivable  that  very  large  transmission  losses 
(2  to  3 dB)  could  be  tolerated  during  a homing  mission  and  not 
degrade  target  intercept.  An  increase  in  radome  transmission 
losses  is  caused  by  three  phenomena: 

1.  A change  in  electrical  thickness  that  produces  a reflec- 
tion loss, 

2.  An  increase  in  the  loss  tangent  of  the  radome  material 
that  produces  an  absorption  loss,  and 

3.  Distortion  of  the  antenna  aperture  illumination  that 
changes  the  antenna  gain  (this  change  may  produce  a gain 
rather  than  a loss) . 

Laboratory  experiments  at  room  temperature  on  a Pyroceram 
C-band  radome  verified  that  reflection  losses  were  less  than  1 uB 
over  a 10%  frequency  bandwidth.  The  ratios  of  the  square  roots  of 
the  dielectric  constants  of  Pyroceram  to  alumina  and  to  fused 
silica  are  0.76  and  1.31,  respectively.  Using  these  ratios,  1-dB 
loss  bandwidths  for  alumina  and  fused  silica  are  estimated  to  be 
7.6%  and  13%,  respectively.  Thus,  reflection  losses  are  less  than 
1 dB  for  relative  radome  design  changes  of  less  than  3.8%,  5%,  and 
6.5%  for  alumina,  Pyroceram,  and  fused  silica,  respectively.  Also, 
it  has  been  found  that  absorption  loss  (in  dB)  of  a half-wave  ra- 
dome (normal  incidence)  for  loss  tangents  (tan  6)  of  less  than  0.5 
can  be  approximated  by  13.65  tan  6.  This  expression  indicates  that 
tan  6 can  be  as  high  as  0.05  before  0.7  dB  of  absorption  loss  is 
realized.  Therefore,  in  this  study  we  may  safely  assume  that  ra- 
dome transmission  losses  will  have  little  influence  on  the  electri- 
cal limitations. 

Angular  Error  Prediction  Method 

An  analytical  method  was  developed  for  predicting  the  ra- 
dome angular  error  slope  (e)  during  flight.  This  method  is  based 
on  well-known  mathematical  relations  and  experimentally  determined 
constants.  The  experimental  and  mathematical  characteristics  as- 
sumed for  the  radome  model  are  listed  below: 

1.  The  radome  has  a von  Karman  shape  with  a length-to- 
diameter  ratio  of  2.1. 
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2.  The  antenna  is  an  11-in. -diameter,  four-quadrant  array 
with  an  effective  phase-center  separation  of  3.89  in. 

3.  The  RF  bandwidth  is  zero;  i.e.,  single-frequency  opera- 
tion. 

4.  The  radome  has  been  corrected  to  have  zero  error  slope 
at  room  temperature. 

5.  The  basic  thickness  of  the  corrected  radome  in  wave- 
lengths is 


" Wecr  8in*6)  * (43) 

where 

8 * tan  ^ £q  • Brewster's  Angle, 

and 

dQ  is  the  room  temperature  radome  wall  thickness, 

Aq  is  the  center  design,  free-space  wavelength  of 
Incident  radar,  and 

eg  is  the  room  temperature  dielectric  constant. 

6.  The  worst-case  average  error  slope  occurs  in  the  gimbal 
region  through  the  radome  nose  and  is  linear  between 
±10°.  This  average  error  slope  is  derived  from  experi- 
mental data  by  averaging  the  worst-case  angular  error 
over  +10* ; i.e. , 


e 

max 

20° 


■ Ae 


7.  The  worst-case  angular  error  slopes  are  associated  with 
polarizations  parallel  to  the  plane  of  rotation;  i.e., 
the  E-plane. 


- 35  - 


THE  JOHNS  HOPKINS  UNIVERSITY 

APPLIED  PHYSICS  LABORATORY 
laurel.  Maryland 


8.  Radome  transmission  losses  resulting  in  lower  receiver 
sensitivity  do  not  degrade  missile  homing  performance. 

9.  A change  in  electrical  thickness  is  directly  related 
to  a change  in  angular  error  slope;  i.e., 


— * C 7^ 
e0  1 *0 


9 


where 


* 


2’ffo  - 8l"  * 


(i|  is  the  electrical  thickness  of  the  radome; 
$ is  the  radar  incidence  angle). 


Laboratory  measurements  at  C,  X,  and  K bands  have  indicated 
that  angular  error  magnitudes  are  approximately  inversely  propor- 
tional to  the  antenna  diameter  in  wavelengths  or  to  interferome- 
ter gain.  Also,  other  laboratory  investigations  with  a low-di- 
electric-constant radome  sandwich  design  have  been  made  (Ref.  13) 
and  the  results  indicate  that  angular  errors  are  generally  lower 
over  a broadened  frequency  band  for  materials  of  lower  dielectric 
constant.  In  particular,  data  obtained  on  a C-band  Pyroceram 
radome  and  a two-ply  sandwich  wall  radome  indicated  that  the  im- 
provement in  frequency  bandwidth  is  approximately  proportions*  to 
the  ratio  of  the  square  roots  of  the  dielectric  constants  of  the 
two  radome  deolgns.  A mathematical  expression  describing  these 
results  is 


Ae  w (Dielectric  Constant  Function)8 
e^  (Interferometer  Gain)  <Jig  * 


Ref.  13. R.  H.  Hallendorff,  'Videband  Radome  Antenna  Re- 
search," Section  3/9,  "Research  and  Development  Quarterly  Report," 
APL/JHU  U-RQR/64-3,  July-September  1964. 
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which  can  be  written  as 


AVEQ  ~ 1 Mi 

4»  D ’ 


or 


Ae  * 


(44a) 


(44b) 


where  D is  the  antenna  diameter  and  B * AegMu.  The  eg  term  in 
the  numerator  of  Eq.  (44a)  is  included  for  converging  the  error 
magnitude  to  zero  for  the  case  of  a radome  material  having  a rela- 
tive dielectric  constant  of  1.  The  constant  B in  Eq.  (44b)  was 
determined  empirically  from  error  slope  data  taken  at  C-,  X-,  and 
K-band  frequencies  on  Pyroceram  radomes.  Experimental  average 
angular  error  slopes  occurring  through  the  radome  nose  were  de- 
termined as  a function  of  frequency  and  plotted  as  a function  of 
percent  change  in  frequency  (Af/fQ).  The  center  frequency  (fg) 
is  defined  as  the  frequency  where  the  error  slope  through  the  nose 
equals  zero. 

Figures  4,  5,  and  6 show  the  typical  change  in  experimental 
average  error  slope  through  the  nose  as  a function  of  percent 
change  in  frequency  for  a constant-wall  half-wave  radome.  The 
K-band  data  of  error  slope  versus  frequency  change  shown  in  Fig.  5 
are  estimated  from  boresight  error  measurements  that  were  made  on 
a full-wave  (two  half-wave  thicknesses)  and  a third-order  (three 
half-wave  thicknesses)  radome  using  a pair  of  standard-goin  horns 
with  a gain  approximately  10  dB  down  from  that  expected  from  an 
11-in. -diameter  antenna.  Also  shown  in  these  figures  are  straight 
line  approximations  to  the  experimental  data.  The  slopes  of  the 
experimental  functions  generated  in  Figs.  4 and  5 are  shown  to  be 
equal  to: 
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Fig.  5 X-and  K-Band  Borasight  Error  Slope  versus  Percent  frequency 
Change 


! 


Boresight  error  slope  change  per 
1%  change  in  frequency 
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Fig.  6 Comparison  of  Error  Slope  Magnitude  versus  Reciprocal  of  Antenna  Aperture 
Sue  in  Wavelengths  (2X0/D)  ''  (Pyroceram  9606) 
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The  numerical  expression  derived  from  experimental  data  for  use 
in  the  radome  analysis  is 


Ae  * 1.615 


deg/deg 


(46) 


where  a ■ D/2.  An  increase  in  error  slope  (Ae)  can  be  determined 
from  Eq.  (46)  for  a given  change  in  electrical  thickness  (Aiji/^o) 
once  the  frequency  (wavelength),  relative  dielectric  constant,  and 
effective  phase  center  separation  (a)  are  known. 

What  remains  now  is  to  calculate  A^/ipo  for  the  radome  wall 
as  a function  of  trajectory  time.  To  do  this  the  wall  is  assumed 
to  be  divided  into  m subsections  through  the  thickness.  Any 
single  section  (n)  is  considered  to  have  its  own  value  of  electri- 
cal thickness  (ijin) : 


* 


n 


2nd 
n 

A 


- sin  $ 


(47) 


and  the  average  electrical  thickness  of  the  whole  wall  is: 


*1* 

avg 


1 

m 


(48) 


The  relative  electrical  thickness  change  is  then 


££  « -£Y£  I _ ! 


♦r 


*p 


n«m 


-E 


n«l 


4 


eA  + Ae  - sin  $ 
0 n _ 

eQ  “ sin  $ 


(*•?) 


(49) 


- 41  - 


TH£  JOHNS  HOPKINS  UNIVERSITY 

APPLIED  PHYSICS  LABORATORY 

uuoei.  Maryland 


In  the  URLIM  program  the  electrical  thickness  change 
(Eq.  (49))  is  calculated  using  the  same  wall  thickness  divisions 
used  for  the  thermal  stress  analysis  in  subroutine  SIGMET.  The 
result  is  made  available  to  routine  MOBSER  where  Eq.  (46)  is 
solved  for  the  boresight  error. 


AERODYNAMIC  LOAD 

During  the  flight  of  a missile,  considerable  aerodynamic 
pressure  and  skin  friction  can  develop  on  the  radome  because  of 
its  high  speed  and  angle  of  attack.  Also,  during  maneuvering  to 
the  target,  high  lateral  and  angular  accelerations  can  be  devel- 
oped, causing  inertial  forces  in  the  radome  structure.  Calculat- 
ing the  magnitude  of  these  forces  ana  the  resultant  stresses  is 
important  in  radome  design  in  order  to  determine  the  limitations 
of  the  radome  and  hence  the  missile.  The  calculation  procedure 
is  designed  to  be  applied  at  all  times  during  the  flight  of  the 
•"iasile.  Moreover,  if  the  calculations  are  done  over  a variety 
of  trajectories  that  cover  the  missile's  propulsion  capabilities, 
then  the  mechanical  load  limits  for  the  radomo  can  be  accurately 
stated. 


The  present  analysis  of  mechanical  loadB  on  a radome  con- 
siders forces  caused  by  axial  and  normal  pressure  drag,  axial 
friction  drag,  normal  (or  lateral)  and  axial  accelerations,  and 
angular  accelerations.  Figure  7 is  a schematic  of  a radome  with 
the  various  forces  listed  that  can  act  at  any  time  during  a flight. 
F{jr  and  F^r  are  the  normal  and  axial  aerodynamic  forces  acting  on 
the  radome  and  are  shown  acting  at  the  radome 's  center  of  pres- 
sure (cp).  The  inertial  forces  are  shown  acting  at  the  center  of 
mass  of  the  radome,  mr.  The  sum  of  these  forces  can  be  resolved 
into  forces  and  moments  acting  at  the  base  of  the  radome  (point  0, 
Fig.  7).  These  forces  in  turn  cause  stresses  in  the  walls  of  the 
radome  that  can  be  calculated  from  knowledge  of  the  radome 's  base 
radius  (RB)  and  wall  thickness  (t). 

Figure  7 is  a schematic  of  a radorpe  ^uring  flight  with  three 
coordinate  systems.  Coordinate  system  X ,Y  is  inertial  and  fixed 
to  the  earth;  system  X,Y  is  centered  at  the  center  of  mass  of  the 
missile  and  is  aligned  at  any  instant  with  the  flight  path  (i.e., 

X is  along  the  velocity  vector  and  Y is  perpendicular  to  it). 
Lateral  maneuvers  are  considered  to  be  along  the  Y-axis.  The 
third  coordinate  system  (x,y)  is  J coated  at  the  radome 's  base  and 
is  fixed  to  the  missile  body;  i.e.,  it  is  aligned  with  the  missile 
axes,  x along  the  centerline  and  y perpendicular  to  the  centerline 
at  0.  All  directions  shown  in  Fig.  7 are  positive. 
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X',Y'-  Coordinate  system  attached  to  ground 

X,Y-  Coordinate  system  along  and  perpendicular 

to  velocity  vector  of  missile  at  center  of  mass  of 
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During  the  flight,  the  aerodynamic  pressure  and  friction 
forces  act  on  the  radorae  with  a magnitude  and  direction  that  is 
dependent  on  speed  (Mach  number)  and  angle  of  attack  (ot) . For 
various  shapes,  these  forces  are  tabulated  in  coefficient  form 
and  are  divided  into  axial  (x-direction)  and  normal  (y-direction) 
components,  as 


and 


Nr 


<VoMo 


Ar 


<VoMo 


/2)A 


/2)A 


(50a) 


(50b) 


where  CAr  and  CNr  are  the  force  coefficients  for  the  radome  in  the 
axial  and  normal  directions;  YqPqMq2/2  is  the  dynamic  pressure  of 
the  stream  with  static  condition  or  pressure  pq,  velocity  (Mach 
number)  Mq,  and  specific  heat  ratio  Y0;  and  A is  the  area  of  the 
base  of  the  radome.  These  component  pressure  forces  are  said  to 
act  at  the  center  of  pressure  (cp)  where  their  effect  on  the  ra- 
dome  is  resolved  into  the  component  forces  with  no  moment.  The 
location  of  this  point  itself  depends  on  Mach  number  and  angle  of 
attack  and  may  be  positioned  away  from  the  body  centerline.  As 
shown  in  Fig.  7,  the  location  of  the  center  of  pressure  is  defined 
by  xCn  and  yCp-  Reference  14  is  an  excellent  source  for  these  co- 
efficients for  several  nosecap  shapes.  For  the  present  analysis, 
the  center  of  pressure  is  assumed  to  vary  only  along  the  x-axis 
and  to  have  no  displacement  from  the  centerline  (i.e.,  ycp  * 0). 

During  an  engagement,  the  missile's  control  system  will  re- 
quire rapid  changes  in  the  flight  path  direction.  To  accomplish 
this,  the  missile  will  change  its  angle  of  attack  to  provide  suf- 
ficient lift  for  the  maneuver.  The  present  analysis  considers 
only  lateral  accelerations  caused  by  angle  of  attack  changes  in 
the  X,Y  plane.  At  some  angle  of  attack  during  a maneuver,  the 
rate  of  change  of  angle  of  attack  may  be  quite  large.  Dynamic 
simulations  of  various  missile  configurations  are  available  to  sup- 
ply these  values.  The  causes  of  inertial  forces  that  can  be  con- 
sidered to  be  acting  at  the  center  of  mass  of  the  radome  are 


Ref.  14.  E.  T.  Mar ley  and  H.  Ginsberg,  "Supersonic  Pressure 
Distribution  and  Axial  Force  Characteristics  of  Axisymmetric  Noses 
at  Angle  of  Attack,"  Paper  Presented  at  the  Seventh  U.S.  Navy  Jym- 
posium  on  Aeroballistics,  7-9  June  1965. 
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• • 


Lateral  acceleration:  mrY, 

•• 


Axial  acceleration: 
Angular  velocity: 


Angular  acceleration:  mrR  \a  + 9 /, 

Radome  mass:  t^g. 

These  forces  are  shown  in  Fig.  1 acting  in 

The  resultant  forces  at  the  base  of 
about  point  0 are  found  to  be 


and 

positive  directions, 
the  radome  and  moment 


F - F.  + m 
x Ar  r 


g sin  (a  + 0)  + X cos  a 


+ Y sin  a - R ( a + 0 )^  ] , 


F - -F„  + m 
y Nr  r 


g cos  (a  + 6)  - X sin  a 

+ Y cos  a + R ( a + 0 ) ] , 

M0  " _FNr  Xcp  + mrxmr  [8  cos  (a  + e)  " * sin  a 
+ Y cos  a + 11  (a  + 0)  j , 


(51a) 


(51b) 


(51c) 


These  stresses  caused  by  the  forces  are  distributed  around 
the  base  of  the  radome  and  are  maximum  either  at  the  top  (y“RB)  or 
at  the  bottom  (y«-RB)  or  the  radome 's  base.  These  stresses  are 


0 

y 


-RB 
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— RB 
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and 


(52a) 


(52b) 
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The  angle  of  the  missile's  velocity  vector  with  the  local 
horizongal  (0)  is  often  called  the  quadrant  elevation.  It  is  de- 
termined from  the  trajectory  information  input  by  the  user,  namely 
the  altitude  and  velocity  as  functions  of  time.  In  the  brief  deri- 
vations that  follow  it  will  be  assumed  that  the  altitude  and  ve- 
locity are  known  at  two  points  along  the  trajectory;  i.e.,  and 
are  known  at  t]_,  and  Z2  and  V2  at  t2*  Further,  it  is  assumed 
that  the  vertical  acceleration  of  the  missile  (Z)  is  constant  dur- 
ing the  time  interval  t2  “ tj  and  that  the  quadrant  elevation  does 
not  change  (i.e.,  0^  ■ 62  * 9)*  If  2 is  constant,  then 


Z - C t + c2  , 

and 

Z - (C1/2)t2  + C2t  + C3  . 

At  the  two  times  t]^  and  t2»  four  relationships  may  be  written; 


Zx  - Vx  sin  6 - C1t1  + C2  , (53a) 

t2  - V2  sin  0 - C1t2  + C2  , (53b) 

Z2  - (C1/2)t22  + C2t2  + c3  , and  (53c) 

Zx  - (C1/2)t12  + C2tx  + C3  . (53d) 

Equation  (53a)  may  be  rewritten  as 

C2  - Vx  sin  0 - , (53al) 


and  substitution  into  Eq.  (53b)  gives 


Vj  sin  0 « C^t2  + sin  0 - . 


(53bl) 
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If  we  assume  ti  » 0,  so  that  t2  will  now  represent  the  time  inter- 
val between  points  1 and  2 then  Eq.  (53d)  becomes 


C 


3 


(53dl) 


Eq.  (53bl)  becomes 


sin  6 - C^t^  + V1  sin  0 , 


(53b2) 


and  Eq.  (53c)  becomes 
C1  2 

Z2  ‘ ? t2  + (V1  8ln  0)t2  + Z1  * 


(53d) 


Equation  (53b2)  may  be  rewritten  ass 

C1  " (V2  “ V(sin  0)/t2  ‘ 

Substitution  of  Eq.  (53b3)  into  Eq.  (53cl)  gives 


e 


2(Z2  - Zx) 

*2*2  +'V 


(53b3) 


(54) 


The  solutions  to  Eqs.  (51a),  (51b),  (51c),  (52a),  and  (52b) 
are  obtained  by  the  use  of  Eq.  (54)  in  the  AERLOAD  subroutine  and 
provide  attachment  stresses  as  a function  of  trajectory  time  in 
the  URLIM  program. 
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2 . NUMERICAL  METHODS 


TRANSIENT  HEAT  CONDUCTION 

Consider  a body  divided  into  N mass  lumps  (nodes)  as 
shown  below: 


In  general,  any  node  n may  be  connected  thermally  to  any  number 
of  other  nodes.  To  determine  the  temperature  of  any  node  (T'  ) 
due  to  Fourier  conduction  within  the  material  at  some  time  in- 
terval away  from  an  initial  starting  time,  the  following  equation 
is  solved  for  T'n: 


pV  c (T' 

- jljlIIjgl 
^r__ 


where 


k A 
n-i  n-i 


n-i 


p is  the  density  of  the  material, 

V is  the  volume  of  the  node, 
n 


(55) 
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c is  the  specific  heat  of  the  material, 

P 

t is  time  (t  being  the  current,  time  and  t being  the 
time  in  the  future), 

k . is  the  thermal  conductivity  between  n and  each 
n”  connecting  node  i, 

A is  the  area  perpendicular  to  heat  flow  between  n 
n and  each  connecting  node  i, 

L is  the  distance  between  n and  the  connecting 

n~  nodes  i,  perpendicular  to  A, 

T is  temperature  (T*  being  at  t*,  T being  at  t),  and 

K is  the  total  number  of  directly  connecting  nodes. 

The  left-hand  side  of  Eq.  (55)  represents  the  change  in  internal 
energy  of  the  node  n during  the  time  interval  t - t.  The 
right-hand  term  of  Eq.  (55)  is  the  net  heat  transferred  to  or 
from  the  node  by  all  the  nodes  it  is  in  contact  with.  Equa- 
tion (55)  is,  moreover,  the  finite  difference  approximation  to 
Eq.  (4a)  derived  earlier  in  this  volume  (Chapter  4 of  Ref.  15 
also  contains  this  derivation).  The  right-hand  side  of  Eq.  (55) 
can  have  additional  terras  which  include  other  modes  of  heat  trans- 
fer between  nodes  as  well  as  the  application  of  boundary  condi- 
tions for  external  surface  nodes.  The  equations  listed  below  de- 
fine how  these  additional  modes  of  heat  transfer  are  calculated 
within  the  program. 

Radiation  Term 

R 

q »o  S % F A e l 7 - T M 

qR  JLd  n-i  n-i  n-i  \ n if 

i-1 

where 

o is  the  Stefan-Boltxman  constant, 

F is  the  view  factor  from  n to  i, 


Ref.  15,  G.  H.  Dusinberre,  Heat-Transfer  Calculations  by 
Finite  Differences,  International  Textbook  Co.,  Scranton,  PA, 
1961. 
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e . is  the  infrared  emissivity  factor  from  n to  i,  and 
n-i 

R is  the  total  number  of  nodes  radiating  with  node  n. 
Constant  Flux  Terms 


- £ ^n-l  Vi 


i-1 


where: 


• // 

q 


n-i 


is  the  independently  specified  rlux  between  n and 


the  independent  source  i per  unit  area,  and 
Q"  is  the  number  of  independent  sources. 
Convective  Flux 


Convective  (aerodynamic)  heating  is  a boundary  condition 
that  is  solved  for  implicity  by  the  SAERO  subroutine  and  therefore 
does  not  enter  into  the  nodal  temperature  solution  being  described 
here.  Convective  flux  boundary  conditions  can  be  specified  by  the 
user  on  as  many  surfaces  of  the  thermal  model  as  are  desired. 

Internal  Generation 


where 

t jh 

q is  the  heat  generated  within  node  n due  to  indepen- 
dent phenomena  (e.g.,  nuclear  fission.  Joule  heating, 
chemical  reaction,  Peltier  effect,  etc.),  which  may 
be  time  dependent. 

Solution  Technique 

The  transient  solution  for  N nodal  temperatures  requires 
choosing  a time  step  (t  - t)  and  solving  the  N equations  of  the 
form  of  Eq.  (55)  for  all  values  of  T_'n.  The  N temperatures  are 
stored  in  the  subscripted  variable  T (dimensioned  from  1 to  N) ; 
each  right-hand  side  of  Eq.  (55)  is  combined  into  the  subscripted 
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variable  Q (dimensioned  from  1 to  N) , and  each  pVncp  is  stored 
in  the  subscripted  variable  C (dimensioned  from  1 to  N) . The 
N equations  are  then  reduced  to: 


T'  - T - (Q/C)  (t'  - t)  , 


(56) 


where  the  bar  indicates  a vector  of  dimension  N. 

For  completeness,  and  to  appreciate  the  nonlinearity  of 
Eq.  (55)  or  (56),  it  is  noted  that  the  coefficients  in  these  equa- 
tions are  dependent  on  the  temperature  of  the  node;  l.e.,  k,  e, 
and  cp  are  treated  as  functions  of  temperature.  If  k ■ f(T)  in 
general,  then  for  Eq.  (55): 


and  if  e ■ g(T)  in  general,  then 


The  special  case  of  two  nod  js  of  differing  materials  (hence  dif- 
ferent k functions)  connected  together  through  a contact  re- 
sistance is  handled  by  the  COMCON  subroutine,  which  evaluates  the 
following  equation: 


cc 


■E 

i-1 


cc 


<Tn  - V 


where 


(59) 


(59a) 
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Subscript  x denotes  the  location  of  the  interface  between 
nodes  n and  i,  and 


^2  * ^2  ’ t*le  t^ierma^  conductance  of  the  contact,  and 


Implied  in  Eq.  (59)  is  a solution  for  Tx  that  is  carried 
out  by  subroutine  COMCON  and  is  based  on  the  fact  that  the  inter- 
face has  no  thermal  capacity  (therefore  all  the  heat  conducted  to 
one  side  of  the  contact  is  completely  conducted  out  at  the  other 
side) . 


The  transient  solution  is  obtained  by  choosing  a time  step 
(t  - t)  and  solving  all  N equations  of  the  form  of  Eq.  (56) 
for  T . One  must  be  careful  in  choosing  the  time  step  so  that  the 
solution  will  be  numerically  stable.  The  rule  for  choosing  a time 
step  is 


At  * 0.9 


MIN 

n*l,N 


(60) 


where  Hn-i  is  the  net  conductance  between  nodes  n and  i and 
can  be  composed  of  conduction  and/or  radiation  terms: 


n-i 


k . A . 
n-1  n-i 

L v 

n-i 


+ Tn)(Tl2  + T»2)  ' 


MIN(  ) represents  a choice  of  the  minimum  value  of  all  N such 
values.  (The  N values  of  are  stored  in  the  vector  H 

within  the  computer  program.)  Equation  (60)  can  be  viewed  as  a 
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statement  of  the  second  law  of  thermodynamics;  that  Is,  for  the 
node  in  the  network  with  the  combination  of  smallest  thermal 
capacity  (Cn)  and  the  largest  total  thermal  conductance  (EHn_i), 
the  amount  of  heat  that  would  flow  out  of  the  node  during  At  if 
all  surrounding  nodes  were  at  0°R  would  not  cause  the  tempera- 
ture of  that  node  to  be  less  than  0°R.  Algebraically,  the  rela- 
tionship can  be  seen  by  letting  and  T n approach  0 in  Eq.  (55). 
In  Eq.  (60)  we  have  arbitrarily  taken  90%  of  this  minimum  to  en- 
sure stability. 

In  the  computer  program,  the  STEP  routine  computes  At  and 
the  nodal  temperatures  via  Eqs.  (60)  and  (56),  respectively. 

Orthogonal  Geometries 

The  fundamental  heat  conduction  equation  was  described  in 
Section  1 of  this  report.  In  those  discussions,  the  spatial  coor- 
dinate system  used  in  the  partial  differential  equations  was  re- 
quired to  be  orthogonal.  When  the  partial  differential  equations 
are  approximated  by  difference  equations  (e.g.,  Eq.  (55))  proper 
care  should  be  exercised  to  see  that  the  nodal  network  (l.e.,  the 
geometry)  is  orthogonal.  In  Eq.  (55),  orthogonality  is  incorpo- 
rated in  the  two  terms  An-i  and  L,^,  where  A is  a "contact" 
area  between  adjacent  nodes,  and  L is  the  distance  perpendicu- 
lar to  A between  node  centers.  In  coordinate  systems  with 
natural  curvature  (cylindrical  and  spherical,  for  example)  the 
simple  A/L  approximation  is  not  correct  because  the  effects  of 
curvature  are  ignored.  In  the  following  derivatives,  the  partial 
differential  equations  will  be  written  and  the  characteristic  solu- 
tions for  steady  state  will  be  obtained  for  two  adjacent  nodes. 

The  appropriate  correction  terms  will  then  be  identified. 

Cylindrical  Nodes 

Considering  heat  conducted  in  the  radial  direction  only  in 
the  steady  state,  the  Fourier  heat  equation  (c.f.  Eq.  (2))  is 
written  as 


32T  1 _3T 

3 2 r 3r 

r 


0 . 


If  p * dT/dr,  then  by  substitution  and  allowing  the  partial  deriva- 
tives to  be  total  derivatives. 


dp/dr  + p/r  » 0 . 
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Further,  if  p is  of  the  form  p - dT/dr  * C/r,  then 


T(r)  * C In  r + D , 

where  C and  D are  constants  to  be  determined  by  the  boundary 
conditions  depicted  in  the  following  sketch: 


In  the  sketch  the  nodal  interface  Is  at  r,  the  inner  node  center 
is  at  r*  with  temperature  T^,  and  the  outer  node  center  is  at  r0 
with  temperature  T0.  Using  these  two  conditions  in  the  above 
equation  for  T(r)  will  yield: 


<T  - T ) 

T(r)  ■ To  + In  (r  /r.)  (ln  r * ln  ro>  ' 
o i 


The  heat  conducted  through  the  nodal  contact  is  given  by 


q(r)  - -A(r)  k VT  , 
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where  V is  the  gradient  operator  (c.f.  Eq.  (6))  and  reduces  to,. 
dT/dr  in  this  unidirectional  case.  A(r) * the  function  relating 
cross-sectional  area  to  radius,  for  the  nodes  shown  is 


A(r)  » (61  - 02)rz  , 

where  z is  the  depth  or  axial  dimension,  and  0i  and  02  are  the 
angular  locations  of  the  node  centers.  When  the  above  relation- 
ships for  A(r)  and  T(r)  are  substituted  in  the  equation  for  q(r) 
the  following  equation  results: 


q(r) 


-z  (02  - 02)  k<To  - Ta) 


1°  <r0/ri) 


Comparing  this  equation  with  Eq.  (55)  will  reveal  that  if 


(0L  - e2)  z 
ln  <ro/rt) 


t 


then  proper  consideration  for  the  cylindrical  curvature  will  have 
been  made.  The  curvature  effect  derives  from  the  cross-sectional 
area  having  a dependence  on  the  coordinate  perpendicular  to  the 
area.  In  cases  where  the  area  is  invariant  in  the  direction  of 
heat  flow  (rectilinear  coordinates,  for  example)  no  "correction" 
term  will  be  required. 

For  conduction  in  the  circumferential  direction  in  a cy- 
linder (i.e.,  in  the  0 direction),  the  following  partial  differ- 
ential equ’*-ion  is  written: 


2 

r 


3 2 T(9) 
2 

3 0* 


0 , 


which  simplifies  to 


dT/d0  - C , 
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or 

T(8)  - C6  + D , 

where  C and  D are  again  constants  of  integration  determined 
by  the  boundary  conditions  sketched  in  the  following  figure: 


The  heat  conducted  circumferentially  between  the  nodes  in  the 
sketch  is  given  by 


4(e) 

with 

A(6) 

hence 

q(0) 


A(e>  kin 


(rQ  - rx)  z (T(62)  - T(61))  k 
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where  r is  the  average  radius,  or  the  radius  of  the  node  centers; 
i.e. , 


r - (r  + r.)/2  • 

0 I 


In  this  case  the  term  (rQ  - r^)z  is  the  cross-sectional 
area  and  r($2  - 8^)  is  the  path  length  between  node  centers  and 
there  is  no  correction  term.  This  is  as  expected  because  there 
is  no  variance  of  the  area  normal  to  the  circumferential  direction. 
By  a similar  argument,  the  axial  cylindrical  direction  will  have 
no  correction  term  associated  with  it. 

Spherical  Nodes 

In  spherical  geometries  the  coordinate  system  consists  of 
(a)  the  co-latltudinal  angle  ($)  measured  from  the  pole  of  the 
sphere  with  positive  radial  values  (i.e.,  the  "North"  pole); 

(2)  the  longitudinal  angle  (6)  measured  as  projections  in  the 
equatorial  plane;  and  (c)  the  radial  distance  (r)  measured  from 
the  center  of  the  sphere.  In  the  nodal  sketches  that  follow, 
these  directions  will  be  illustrated  for  clarity. 

Considering  only  conduction  in  the  co-latitudinal  direction 
at  steady  state,  the  Fourier  equation  is 


d2T/d$2  + ctn  $ dT/d*  - 0 . 

If  we  let  dT/d$  « p ■ G esc  then 


dp/d$  + ctn  $ p ■ 0 


(61) 


by  substitution.  Also, 


dp/d$  » -C  esc  $ ctn  $ , 


so  that  by  substitution  in  Eq.  (61), 


-esc  $ ctn  $ + ctn  $ esc  $ ■ 0 , 
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which  is  an  identity  and  verifies  the  choice  of  p - C esc  4> . 
Hence, 


dT/d<j>  * C esc  $ , 


and 


T<$)  - C In  (tan  $/2)  + D , 


where  C and  D are  integration  constants  determined  by  the 
boundary  conditions  described  next.  The  sketch  below  shows  two 
spherical  nodes  emphasizing  the  co-latitudinal  coordinate  (<j>) : 


In  the  following,  » T(4>^),  T2  » TC^),  and  the  node  centers  are 
located  at  ($1,  r,  0)  and  (d>2 » r»  0)  * where  0 is  the  average 
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meridional  coordinate  ((0|  + 02)/2).  The  area  normal  to  the  4> 
direction  is  given  by 


AU) 


sin  d>  (rQ2  - r12)  (02 
2 


» 


where  r0  is  the  radius  of  the  outer  surface  of  the  nodes  and  r^  is 
the  radius  at  the  inner  surface  of  the  nodes;  hence  r0  - r^  is 
the  radial  thickness  of  the  nodes.  Using  only  the  $ direction 
terras  for  the  spherical  gradient  operator  (c.f.  Eq.  (7)), 


q($)  - -A($)  k (1/r)  dT/d*. 


The  previous  relationship  for  T($)  combined  with  the  boundary  con- 
ditions will  give 

dI  .tn  ♦ ;t2  - Tt) 
d4  (cos  <^2  ~ coa  ♦j^)  * 


Combining  the  A($)  and  dT/d$  relationships  will  yield 

r±)  (02  - 6^  k (T2  - Tx) 

'tan  (♦2/2)]  * 

ln  tan  (♦1/2) 

giving  an  effective  A/L  in  the  $ direction  of 


A 

L 


* 


- V (°2  - V 

tan  (<p2/2)  1 
tan  (4^/2) 


When  considering  only  heat  conducted  in  the  radial  direc- 
tion, the  divergence  of  the  temperature  field  (T)  for  spherical 
geometries  can  be  reduced  to 
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d2T/dr2  + (2/r)  dT/dr  - 0 . 


To  find  T(r).  let  p * dT/dr,  so  that 


dp/dr  + 2p/r  ■ 0 . 

2 

Further,  if  p ■ C/r  , then  by  substitution  in  the  previous  equa- 
tion 

-2Cr3  + (2/r)  (C/r2)  - 0 , 


which  is  an  identity,  verifying  the  choice  p * C/r2  » dT/dr.  One 
integration  of  this  relationship  gives 


T(r)  - (-C/r)  + D , 


where  C and  D are  constants  determined  by  the  boundary  condi- 
tions imposed  as  illustrated  in  the  following  sketch: 
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In  the  sketch,  two  spherical  nodes  are  in  contact  at  radius  r 
and  have  node  centers  at  r0  and  r^.  The  temperature  (which  de- 
pends only  on  the  dimension  r)  is  T0  at  r0  and  at  r^.  Sub- 
stituting these  values  as  boundary  conditions  gives 


and 


dT 

dr 


The  area  normal  to  the  radius  direction  is  given  by: 


2 

A(r)  » r (cos  ^ - cos  $2>  (62  - 0^  . 


The  heat  conducted  across  the  area  at  radius  r is  then 


$(r)  - -A(r)  dT/dr 


(cos  ♦1  - cos  $2)  (62  - 0^)  k (T2  - T^) 
1_  _ 1_ 


By  inspection  with  the  rectangular  case,  the  effective  A/L  for 
this  condition  is: 


A 

L 


r 


(cos  4^  - cos  $2)  (e2  - e ) 
1 1 


When  the  only  dependence  of  temperature  is  in  the  longitu- 
dinal direction  (6)  then  the  divergence  operator  on  T yields 
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CSC  $ 


0 , 


which  reduces  to 

d2T/d92  - 0 , 


and  is  identical  to  the  rectangular  case.  Moreover,  it  is  apparent 
that  the  area  perpendicular  to  the  0 coordinate  does  not  vary  with 
8 as  can  be  seen  in  the  attendant  sketch: 


From  the  simple  form  of  the  divergence  of  T (stated  above)  , by 
inspection 


dT/d0 
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with  being  the  nodal  temperature  at  and  To  being  the  nodal 
temperature  at  02*  The  contact  area  between  the  nodes  is  at  0; 
the  area  perpendicular  co  the  6 coordinate  is: 


A(9) 


(♦2  * V (r„2  - r/> 


y 


and  the  heat  conducted  through  this  area  is  then 


q(9) 


-A  (6) 


k dT 
r sin  $ d0 


The  term  in  brackets  is  seen  as  the  actual  A/L  for  the  meridional 
direction  if  r * (r0  + r^)/2,  ♦ • ($1  + $2^/2,  A is  given  by 
A(6)  above,  and  L is  the  arc  length  (r  sin  $)  (©2  “ ®l). 

In  summary,  by  taking  proper  consideration  for  curvature 
in  cylindrical  and  spherical  coordinate  systems,  effective  A/L 
formulas  have  been  derived  that  should  be  used  in  the  finite  dif- 
ference equations  for  heat  conduction.  The  terms  derived  reduce 
to  the  correct  partial  differential  equations  in  the  limit  of  in- 
finitely small  differences. 


STEADY-STATE  TEMPERATURE  FIELDS 
Iterative  Solution 


In  the  steady  state  there  is  no  change  of  temperature  from 
one  time  to  the  next  so  that  T • T in  Eq.  (56)  and  we  have  no 
net  flow  of  heat  at  any  given  node;  i.e,, 


Q - 0 , 


(62) 


which  states  that  the  net  heat  flow  for  each  node  must  be  zero. 
Solving  the  steady-state  problem  then  involves  finding  the  set  of 
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temperatures  (T)  so  that  Eq.  (62)  is  satisfied.  Since  all  of  the 
individual  terms  that  comprise  Eq.  (62)  are  temperature  dependent 
(c.f.  Eqs.  (57)  through  (59)),  then  Q is  a function  of  T;  i.e., 


Q - F(T)  . 


(63) 


Newton's  method  for  finding  the  zeros  of  a function  can  be 
used  to  solve  this  system  of  equations  in  the  following  way:  If 

an  initial  guess  for  T,  say  To,  is  substituted  in  Eq.  (63),  then 

Qq  “ F(T0)  i 0 . (64) 


A set  of  correction  terms  (Sg)  will  exist  so  that  Tg  - 6g  ■ T and 


F(T)  - 0 , 
or 

f(Tq  - 6q)  - 0 . (65) 


Equation  (65)  is  expanded  in  a Taylor's  series  as 


- _ , _ «g2  F'(T0> 

F(T0)  - 60F  (T0)  + — 21  “ 


0 , 


(66) 


where  the  primes  indicate  derivatives  with  respect  to  temperature. 

It  is  assumed  that  the  initial  guess  (Tg)  is  sufficiently 
accurate  so  that  the  corrections  (6g)  are  small.  Further,  the 
second  and  all  higher  derivative  terms  are  assumed  small  compared 
with  the  fijrst  derivative  term.  With  this  assumption  Eq.  (66)  is 
solved  for  fig: 


7 F(V 
» _ ■ 

° F< (T0) 


(67) 
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Because  of  the  approximations,  Eq.  (67)  can  be  written  as 


and  moreover, 


(68) 


(69) 


where  6^  is  an  approximation  to  the  true  correction  terms  (6q) , 
and  is  an  approximation  to  the  true  steady-state  temperatures 
T.  Combining  Eqs.  (68)  and  (69)  yields: 

T1  - TQ  - [F(T0)/F'(T0)) 


An  iterative  process  can  be  devised  wherein  each  set  of 
temperatures  are  substituted  back  into  Eq.  (62),  new  values  for 
F(Tj^)  and  F (T^)  are  used  in  Eq.  (70)  to  find  a further  improved 
set  of  temperatures  (T2).  The  process  is  repeated  in  this  manner 

until  the  corrections  (6^)  approach  zero  and  the_steady-state  heat 
f lows  (Q)  equal  zero.  Error  tolerances  for  the  6^  values  on  the 
order  of  ±1*R  nave  been  used  successfully  in  practice;  the  precise 
value  is  selectable  when  the  program  Is  executed. 

To  accomplish  the  steady-state  solution  in  this  manner,  it 
is  required  to  have  an  initial  guess  (Tq)  and  to  evaluate  the  _N 
correction  factors  (Qo/Q'o^*  The  heating  rate  for  each  node  (Q) 
is  calculated  as  discussed  above  for  the  transient  solution  and 
the  derivatives  (Q1)  are  calculated  as  follows: 
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Conduction 


ST 


3 

ST 
n n 


£ r 


i-1 


n-i 

K 


(T  - T.) 
n i 


EtN  /£ 

i-1  n-i  i-1 


(T  - T ) — 
UN  i'  L 


Sk 


n-i 


ST 

n-i  n 


(71) 


Radiation 


°_\m  _i. 

ST  0 ST 
n n 


V F , A , e . (t4  -t.M 
n-i  n-i  n-i  \ n i / 


i-1 


Y'  F A e 4T  - 
/ j n-i  n-i  n-i  n 


i-1 


+ V'  F A (t  4 - T M 3cn-i 
l-J  n-i  n-i  \ n if  aT 

i-1  n 


(72) 


F.xtemal  Fluxes  and  Internal  Heat  Generation 


3 £ 'V. 

i-1 

ST 

n 


and 


Sc 

ST 


n 


(73) 


by  definition. 
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Composite  Conduction 


3q 

cc 

3T 


K (T  - T.) 
cc  n i 


i-1 


3K 


3T 


(T  - Tj 


cc 


L*  V 3T 

i-1 


E '< 


i-1 


. (74) 


It  can  be  seen  from  Eq.  (57)  that 


3k  f\  (T  + T )/2 
n-1  In  i 

3T  2 


(75) 


and  similarly,  from  Eq.  (58), 


3e„  4 (T  + T,)/2 

n-1  l n l 

3T  " 2 


(76) 


With  Eqs.  (59a),  (75),  and  (76)  substituted  into  Eqs.  (71) 
through  (74),  the  iterative  process  implied  by  Eq.  (70)  may  be 
carried  out  to  arrive  at  the  steady-state  temperature  distribution. 

Usage  experience  with  this  iterative  solution  technique  has 
been  varied.  In  general  the  following  cautions  are  offered 
(1)  The  number  of  iterations  required  for  solution  depends  strongly 
on  the  initial  guess  for  the  temperature  field  (Tq);  (2)  the 
larger  the  number  of  nodes,  the  more  iterations  will  be  required 
and  the  smaller  the  probability  of  satisfactory  convergence;  and 
(3)  a constant  relaxation  factor  (described  below)  may  improve  con- 
vergence speed  and  solution  accuracy. 

Relaxation  Factor 


A provision  has  been  added  that  allows  Eq.  (70)  to  be 
written  as 


T 


1 


(77) 
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where  R is  a constant  relaxation  factor.  Simply  stated,  values 
of  R greater  than  1 tend  to  accelerate  the  convergence  to 
solutions  that  otherwise  require  many  iterations  and  values  of  R 
less  than  1 tend  to  damp  succei  aive  solutions  to  problems  that 
tend  to  oscillate.  Oscillation  cau  be  characterized  by  successive 
values  of  6 being  positive  and  then  negative,  and  not  diminishing 
in  absolute  value.  The  proper  choice  for  R is  a subject  of  ex- 
tensive mathematical  significance  and  is  generally  found  by  experi- 
ence (i.e,,  trial  and  error). 

The  foregoing  comments  serve  to  Indicate  an  Inherent  unde- 
sirability in  the  iterative  technique  just  described  — its  unpre- 
dictable behavior.  To  avoid  this,  another  scheme  is  available  — 
the  implicit  technique  — and  is  described  next. 


Implicit  Technique 


The  right  side  of  Eq.  (55)  describes  heat  transferred  via 
Internal  conduction,  but  there  also  may  be  radiation  transfer,  as 
well  as  independently  specified  external  heat  fluxes  or  internal 
heat  generation.  Terms  such  as  these  enter  Eq.  (55)  either  as 
constants  or  terms  that  may  have  some  dependence  on  Tn.  The  spe- 
cific form  of  these  terms  is  described  above  for  the  transient 
solutions.  It  will  be  sufficient  for  purposes  of  description  to 
examine  the  internal  conduction  term  as  written  in  Eq.  (55)  and 
note  that  che  other  terms  mentioned  simply  add  to  the  coefficients 
in  the  equations  that  will  be  developed. 


Note  that  in  the  steady  state,  T*n  equals  Tn  (there  is  no 
change  in  .temperature  with  time)  and  the  left  side  of  Eq.  (55)  is 
zero.  Equation  (55)  is  expanded  to 


E 


k , A , 
n-i  n-i 


n-i 


T 

n 


k . A . 
n-i  n-i 


n-i 


0 . 


(78) 


For  every  node  in  the  network  an  equation  of  the  form  of  Eq.  (78) 
can  be  written.  For  boundary  nodes  whose  temperatures  (Tj,)  are 
known,  Eq.  (78)  is  not  relevant  and  the  following  form  is  used: 


T^  « constant 


(79) 
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For  a constant  flux  or  Internal  heat  generation  boundary  condition, 
the  heat  flux  in  simply  added  to  the  left  side  of  Eq.  (78)  as  a 
constant.  It  should  be  noted  that  in  the  case  of  convective  heat 
fluxes  a separate  subroutine  (SAERO)  uses  an  implicit  technique  to 
solve  for  the  surface  temperature.  Therefore,  as  far  as  the  nodal 
network  solutions  are  concerned,  convective  surface  temperatures 
enter  as  equations  such  as  Eq.  (79).  It  is  now  clear  that  for  a 
network  of  N node  points,  N equations  of  the  form  of  Eq.  (78) 
can  be  written.  To  find  all  N unknown  temperatures  requires  the 
simultaneous  solution  of  this  system  of  equations.  Written  in 
matrix  form  the  system  of  equations  is 


11 

C12 

C13  ' 

• 

• 

1 

Cln  • 

• 

• 

• CuT 

21 

C22 

C23  ' 

• 

• 

• 

C2n  * 

• 

• 

•.  C2» 

31 

C32 

C33  ' 

• 

• 

• 

C3n  ‘ 

• 

• 

C3K 

o I C 


. X . = 


• CnH  j I T„ 


b 0 0 0 0 0 0 0 0 0 1 0 0 Oil  Tb  1 ConiUat 


Sw J I T*  I I ^fxt 


The  above  system  shows  nodes  1 through  N as  regular  inter* 
nal  nodes,  a boundary  node  with  value  Tj,,  and  node  N with  an  ex- 
ternally  supplied  heat  flux  (Qext^* 

It  should  be  noted  that  the  off-diagonal  coefficients  are 
symmetric;  l.e., 


Ci.J  "CJ.i 


kii ih 
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and  that  the  diagonal  coefficients  are 


k A 

En-i  n-l 
L , 

n—i 

i-i  n 


The  solution  technique  used  to  solve  these  simultaneous 
equations  is  Gaussian  elimination  with  back  substitution.  In  this 
technique  the  equations  are  algebraically  manipulated  so  that  the 
coefficient  matrix  is  reduced  to  upper  triangular  (i.e.,  l's  oc- 
curring along  the  diagonal  with  all  2eros  below  the  diagonal). 

This  is  accomplished,  tor  example,  by:  (1)  dividing  the  first 

equation  by  C ; (2)  multiplying  the  new  equation  1 by  -C21  and 
adding  the  result  to  equation  2,  thus  producing  a zero  for  C21 > 

(3)  dividing  the  new  equation  2 through  by  the  new  C22*  Proceed- 
ing in  this  manner  through  all  N equations  produces  the  upper 
triangular  coefficient  matrix.  Note  that  at  the  end  of  this  pro- 
cedure, the  last  equation  (N)  will  read 


Eq.  No. 


• 

* • 

• 

• 

• 

* 

• 

1 

• 

• 

X 

• 

m 

• 

N-l 

0 

1 

CN-1,N 

T 

N-l 

• 

N 

0 

0 

1 

T 

_ N 

value 

«■« 

which  says  Tj$  - value.  This  result  can  then  be  substituted  into 
equation  number  N-l,  and  T^_^  can  b?  solved  for.  This  back  sub- 
stitution procedure  can  be  followed  up  through  the  reduced  equa- 
tions until  all  N temperatures  are  known. 

While  this  technique  requires  the  storage  of  the  N x N co- 
efficient matrix,  usage  has  shown  that  solutions  are  achieved  much 
quicker  and  are  generally  more  accurate  than  the  iterative  scheme 
described  earlier. 

It  should  be  mentioned  ‘■hat  the  implicit  technique  is  also 
used  iteratively;  i.e.,  successive  solutions  to  the  temperature 
field  are  compared  and  convergence  is  judged  according  to  the  suc- 
cessive differences.  Since  the  values  in  the  coefficient  matrices 
are  temperature  dependent,  successive  coefficient  matrices  will  not 
in  general  be  identical  until  the  temperatures  themselves  do  not 
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change.  Generally,  the  number  of  iterations  required  is  small 
(less  than  10). 


INTERPOLATION  TECHNIQUES 

In  arriving  at  the  transient  or  steady-state  temperature  re- 
sponse of  a thermal  network  the  routines  in  URLIM  make  extensive 
use  of  numerical  interpolation.  Because  all  of  the  pertinent  mate- 
rial properties  are  allowed  to  have  temperature  (and  sometimes  pres- 
sure) dependence,  these  functional  relationships  must  be  modeled. 

The  technique  employed  here  is  the  most  simple  and  direct;  the 
functions  are  represented  by  tables  of  independent  and  dependent 
variable  values  and  linear  interpolation  is  most  frequently  used. 

The  table  look-up  procedure  most  often  employed  is  a simple  serial 
search.  In  general,  the  tables  can  be  as  long  as  required  to  ac- 
curately define  the  particular  function.  Search  time  will  in- 
crease with  increases  in  tabular  length.  The  following  discussion 
will  outline  the  techniques  used  by  the  URLIM  interpolators.  The 
PIF1,  PIF1D,  and  DECIDE  routines  in  particular  use  the  following 
scheme . 


Let  the  independent  variable  be  represented  by  the  vector 
of  N values  X and  the  dependent  variable  by  the  vector  of  N 
values  Y with  the  functional  relationship  between  the  values 
expressed  as. 

Y - f(X)  . 


Each  value  in  Y,  Y*  has  a corresponding  value  in  X,  X^  so  that 


- f(xi)  , 


and  each  pair  X^,  Y^  is  a point  on  the  curve 


y * f(x)  . 


If  a value  of  x is  giver,  (xq)  and  the  corresponding  value  of 
X (yo)  is  required,  then  a search  is  performed  through  the  values 
X until  two  values  are  found: 
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x,  s x„  < x 


i+1 


Note  that  the  vector  X must  contain  values  in  monotonic  increas- 
ing order.  When  the  two  points  that  satisfy  the  inequality  above 
are  found  the  value  yg  is  computed  as 


In  the  cases  where  Xq  < or  xq  > Xjj  (where  Xf$  is  the  last  tabu- 
lar value  in  the  vector  X)  linear  extrapolation  is  employed  using 
the  two  given  end  point  values. 

The  routine  LINLCG  uses  a modified  version  of  the  same 
method  described  above.  In  this  case,  after  the  two  tabular 
values  X^  and  are  found  the  corresponding  values  Y^  and  Y^+^ 

are  replaced  in  the  calculations  by  In  (Y^)  and  In  (Yi+i),  re- 
spectively. The  returned  value  yg  is  then 


exp 


, ..  J1"'1 

ln  Yn 


(x0  - Xt) 


This  result  is  then  a semilogarithmic  linear  interpolation.  For 
values  of  y that  are  negative,  the  sign  is  retained  and  the  ln 
function  is  performed  on  the  absolute  values  of  X.  In  the  spe- 
cial case  of  the  value  of  either  Y^  or  Y^+^  lying  between  +1  and 
-1,  the  logarithmic  functions  are  suspended  and  simple  linear  in- 
terpolation is  performed. 

In  the  BIVLID  routine,  the  functional  dependence  of  a 
variable  on  two  independent  variables  is  modeled.  In  this  case 
the  two  independent  variables  x^  and  X2  have  monotonically  increas- 
ing values  stored  in  the  vectors  Xa  and  Xb  with  individual  values 
Xai  » Xa  and  Xb^  c Xb<  The  dependent  variable  y is  defined  as 

y = f(xa,xb)  and  is  represented  in  the  two-dimensional  ajrray  Y. 

The  subscripts  i and  j are  used  to  index  the  array  Y (a.g., 
Y^j  * f(Xa,Xb)).  The  subscript  i is  used  to  index  the  vector 

Xa  and  the  first  index  position  of  the  array  Y;  similarly,  the 
subscript  j _is  used  to  index  the  vector  X5  and  the  second  index 
position  of  Y. 
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_ When  values  of  xa  and  are  given,  BIVLID  scans  the  vec- 
tors Xa  and  X^,  respectively,  until  the  following  conditions  are 
satisfied: 


(1) 


X S x < X 
®i  3 ai+l 


» 


(2)  X.  S 

bJ 


-S'*!. 


j+1 


The  above  conditions  define  i and  j uniquely.  In  the  case 

where  xa  < Xa. , i is  set  to  1.  Also  when  xa  > , where  n 

i _ n 

is  the  extent  of  the  vector  Xa»  then  i Is  set  to  n.  Exactly 

analogous  conditions  hold  for  the  vector  > ^ and  j . 

With  the  values  i and  j established,  two  intermediate 
values  of  the  dependent  variable  are  found: 


The  required  value  of  y is  then  computed  as 


The  BIVLID  routine  also  returns  either  of  the  partial  de- 
rivatives of  the  function  y depending  on  the  value  of  a control 
code.  In  the  example  shown  above  the  derivative  would  be 
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Similar  expressions  to  the  ones  above  can  be  written  for  cases 
where  3y/3x2  are  required. 

The  routine  has  a facility  for  returning  a value  for  one 
of  the  independent  variables  when  values  for  the  dependent  vari- 
able and  the  other  independent  variable  are  supplied.  This  func- 
tion is  termed  an  "inverse"  interpolation  and  is  useful  in  cer- 
tain applications.  In  this  process  if  the  value  xa  is  given  along 
with  the  value  y(xa,x  fc),  then  the  value  x {,  is  returned.  The 
vector  Xa  is  searched  until 


s x < X 


i+1 


9 


thereby  establishing  the  value  of  i and  two  vectors  of  dependent 
values,  namely  Y.  * and  Yj+j^*.  A series  of  intermediate  depen- 
dent variable  values  are  then  found  between  these  two  vectors; 
i.e. , 


such  that  either 

(1)  Y#x  £ y(x1,x/2)  < Y'2  , 

or 

(2)  y(x1>x<2)  > y'2  . 


- 75  - 


THE  JOHNS  HOPKINS  UNIVERSITY 

APPLIED  PHYSICS  LABORATORY 

laurel  Maryland 


When  either  of  these  conditions  is  true,  then  a value  for  the  in 
dex  j will  be  implicitly  known.  The  returned  valae  of  the  de- 
pendent variable  (x  2)  is  then  given  as: 


A similar  procedure  can  be  outlined  whtn  the  values  y(x 
and  X2  are  given  and  x ^ is  desired. 

If  the  conditions  (1)  and  (2)  are  not  satisfied  after  all 

values  of  j are  used,  then  the  dependent  variable  value  that 

was  given  lies  beyond  the  tabulated  values.  In  such  cases,  the 

routine  will  generate  an  appropriate  warning  message  and  return 

the  last  tabulated  value  of  the  required  independent  variable 

(i.e.,  X„  or  Xb  » where  m (or  n)  is  the  extent  of  the  vector 
“m  n 

Xa  (or  Xb) ) . 

In  cases  where  the  dependent  variable  is  double  valued 
(i.e.,  there  is  more  than  a single  value  of  j that  satisfies 
the  conditions  of  (1)  or  (2)  above),  then  the  smallest  such  value 
of  j will  be  used. 

When  the  data  being  used  for  interpolation  have  an  inher- 
ent logarithmic  dependence,  a linear-logarithmic  interpolation 
may  be  obtained  in  much  the  same  way  as  done  for  the  single  inde- 
pendent variable  data  in  LINLOG,  discussed  earlier.  For  these 
cases  the  BIVLLID  routine  is  used,  which  has  the  same  basic  algo- 
rithm as  discussed  above  for  BIVLID.  The  logarithmic  aspect  is 
accomplished  by  using  the  logarithm  of  the  dependent  variable 
values  rather  than  the  values  themselves.  In  other  words,  where 
Yj^j  is  used  in  the  relations  above,  In  (Y^j)  is  substituted  and 
the  returned  values  are  the  antilogs  of  the  resultant  values. 


INTEGRATION 

Within  the  thermal-stress  routines  (SIGMA  and  SIGMET)  nu- 
merical integrations  are  required  of  various  tabular  functions. 
The  process  is  accomplished  by  the  TRAP  routine  and  essentially 
evaluates  the  Integral  I,  which  is  defined  as 
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I 


F(x)  dx 


The  function  F(x)  is  assumed  represented  in  tabular  form  (i.e.,  a 
vector  of  independent  variable  values  (X^  and  a corresponding  vec- 
tor of  dependent  variable  values  (Y  “ F(X)). 

The  algorithm  divides  the  interval  [a,b]  in  half  and  the 
three  functional  values  y^  » f(a),  y2  * f((a+b)/2],  and  y3  » f(b) 
are  found  by  interpolation.  The  integral  I is  now  approximated 
by  the  area  under  the  three-point  function  via  the  trapezoidal 
rule;  i.e., 


I2  “ ^yi  + y2^  + (y2  + y3)  (b  - a)/4  . 

The  interval  [a,b]  is  now  divided  into  three  parts  and  the  values 

yx  “ f (a) , 
y 2 « f (2a  + b)/3, 
y*  2 " ^(2b  + a)/3,  and 
y4  - f(b) 


are  determined  by  interpolation.  The  integral  is  approximated 
again  by  the  trapezoidal  process  and  is  called  I3.  At  this  point 
the  value  E * 1 - II2/I3I  is  compared  with  a preset  tolerance  value 
of  0.01.  If  E > 0.01  then  the  interval  is  divided  into  fourths,  the 
approximation  I4  is  generated  via  the  trapezoidal  rule,  and  a new 
value  of  E ■ 1 - II3/I4I  is  compared.  This  process  is  continued 
until  either  E £ 0.01  or  20  such  iterations  have  occurred.  If, 
after  20  iterations,  E is  greater  than  0.01  then  an  appropriate 
message  is  printed,  the  ratio  120^19  is  printed,  and  the  value 
I2q  is  used  as  the  value  of  the  integral. 


SIMULTANEOUS  AND  IMPLICIT  EQUATIONS 

Within  the  aerodynamic  heating  routines,  there  is  a require- 
ment for  the  solution  to  implicit  equations  (i.e.  , equations  with 
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terms  or  functions  of  the  unknown  variable  on  both  sides  of  the 
equal  sign  that  cannot  be  solved  in  closed  form).  Also,  in  the 
iterative  steady-state  solution  procedure,  a number  of  implicit 
equations  must  be  solved  simultaneously.  The  numerical  technique 
used  in  both  cases  is  Newton's  method  wherein  successive  correc- 
tions are  calculated  for  the  value  of  the  unknown  quantity  based 
on  partial  derivatives.  A complete  description  of  this  algorithm 
is  presented  in  the  paragraph  entitled  "Steady-State  Temperature 
Fields"  (Section  2) . 

The  set  of  simultaneous  equations  generated  in  the  implicit 
steady-state  method  are  solved  by  Gaussian  elimination  to  achieve 
an  upper  triangular  coefficient  matrix  and  the  unknowns  (tempera- 
ture) are  solved  for  by  back-substitution.  Again,  this  method  is 
described  in  Section  2 and  will  not  be  repeated  here. 
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3.  DATA  MANAGEMENT 


EXTERNAL  STORAGE 

In  the  URLIM  and  SHTP  subroutine  library  it  has  been  found 
useful  to  have  a number  of  variables  that  can  be  accessed  freely 
among  the  subroutines.  This  class  of  variables  is  called  EXTERNAL 
in  the  PL/I  language.  Table  2 of  Volume  2 of  this  report  lists  all 
of  the  EXTERNAL  variables  used  with  brief  descriptions.  Since 
EXTERNAL  variables  must  have  unique  names,  the  list  will  serve  as 
a guide  to  future  changes  in  the  program  that  may  require  addi- 
tional EXTERNAL  variables. 

DYNAMIC  STORAGE 

Many  of  the  arrays  of  data  stored  by  the  URLIM  routines 
are  of  flexible  extent;  i.e.,  their  dimensions  are  program  vari- 
ables and  are  changed  as  required  from  run  to  run.  Many  examples 
are  evident  in  Table  2 (e.g.,  the  nodal  temperatures  (T)  and  the 
steady-state  coefficient  matrix  (HMAT) . As  might  be  expected, 
the  single  most  important  value  that  determines  the  extent  of  the 
various  arrays  is  the  number  of  thermal  nodes  (i.e.,  the  value 
of  variable  LASCAP) . Subroutine  STORE  is  the  routine  that  allo- 
cates the  bulk  of  the  required  dynamic  storage  according  to  the 
parameters  passed  to  it.  In  establishing  the  storage  limits  that 
will  be  required,  one  of  the  important  aspects  to  manage  is  the 
use  by  the  READCP  routine  of  ’'node"  numbers  that  are  beyond  the 
value  LASCAP.  In  describing  the  network  data  it  is  possible  to 
use  an  extended  input  technique  that  will  require  additional  stor- 
age (c.f.  Appendix  F of  Vol.  1 of  this  report).  When  employing 
this  method,  adequate  storage  must  be  allocated  by  the  STORE  rou- 
tine through  the  parameter  XCAPLIM.  Further  explanations  can  be 
derived  from  the  discussion  of  the  STORE  routine  (Appendix  H of 
Vol.  1 of  this  report). 

The  U3e  by  the  SIGMA,  SIGMET,  and  MOBSER  routines  of  nodal 
temperature  positions  in  the  vector  T must  also  be  allowed  for 
if  the  plotting  flag  is  set  in  the  calling  sequence  to  any  of  these 
routines.  This  feature  does  not  require  storage  additional  to  that 
set  by  LASCAP,  but  it  is  imperative  that  the  "node"  number  used  for 
the  thermal  stress  values  not  be  the  number  of  a node  in  the  network. 

Another  important  data  management  technique  employed  within 
URLIM  is  the  use  of  POINTER  variables  to  locate  the  storage  of  the 
various  material  property  tables  and  the  various  time-dependent 
tables.  In  the  case  of  the  thermal  properties,  the  READKK  routine 
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stores  the  memory  locations  of  the  property  data  it  reads  into 
the  EXTERNAL  POINTER  variable  PP.  The  values  in  PP  are  then  made 
available  to  the  routines  that  require  the  thermal  property  tables. 
Indexing  with  the  vector  PP  is  accomplished  via  the  ID  values  sup- 
plied to  the  READRK  routine  (i.e.,  the  material  thermal  property 
codes) . 


For  the  time-dependent  variables  read  in  by  the  READTM 
routine,  the  EXTERNAL  POINTER  array  TMPNTR  is  used  to  store  the 
locations  of  the  tabular  data.  This  storage  information  is  then 
made  available  to  routines  in  the  URLIM  library  for  processing. 

The  indexing  of  TMPNTR,  a two-dimentional  array,  is  first  by  the 
ID  number  used  when  read  by  READTM  and  then  by  the  dependent 
variable  position  in  the  calling  sequence.  To  illustrate,  con- 
sider the  following  hypothetical  call  to  routine  READTM: 

CALL  READTM  (ID#,  TM,  DEP1,  DEP2,  DEP3,  DEP4 , DEP5, 

# ENT,  INFILE): 

where  the  vector  TM  is  the  list  of  time  values  and  the  vectors 
DEP1  through  DEPS  are  corresponding  lists  of  dependent  variable 
values.  Upon  return  from  this  call  to  READTM  the  storage  loca- 
tions of  the  variable  TM  and  DEP1  through  DEP4  will  be  recorded 
in  the  values  of  TMPNTR  (ID#,  1 through  5),  respectively. 

The  Initial  Storage  Area 

One  of  the  "optimising'1  features  of  the  PL/I  optimising 
compiler  is  its  improved  ability  to  mauage  the  dynamic  alloca- 
tion of  storage,  according  to  the  requirements  of  the  particular 
program.  The  computer  code  that  supervises  the  allocation  of 
storage  segments  for  a particular  program  is  supplied  automati- 
cally by  the  optimise?  compiler;  the  only  item  supplied  by  the 
user  is  the  specification  of  the  amount  of  storage  that  will  be 
needed  for  the  dynamic  storage.  This  area  is  termed  the  initial 
storage  area  (ISA)  and  has  a default  value  in  the  present  com- 
piler implementation  of  8000  bytes.  For  most  programs,  and  cer- 
tainly for  the  URLIM  program,  this  is  an  insufficient  size  for 
the  ISA.  Determining  the  proper  ISA  size  is  done  with  the  aid 
of  the  program-generated  storage  report  that  gives  an  accounting 
of  the  actual  storage  requirements  and  the  number  of  times  stor- 
age outside  of  the  ISA  was  required.  Figure  8 is  a reproduction 
of  a typical  storage  report,  the  salient  features  of  which  are 
discussed  in  the  following  paragraph. 

In  Fig.  8 the  size  specified  for  the  ISA  is  given  in  bytes; 
the  amount  of  actual  PL/I  storage  needed  by  the  job  is  given  and 
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the  amount  of  storage  used  outside  o.  the  ISA  is  given.  Further, 
in  cases  where  storage  outside  the  ISA  is  required,  the  number  of 
requests  to  the  system  (OS)  for  such  storage  is  indicated.  In  a 
program  with  a properly  specified  ISASIZE,  the  amount  of  storage 
required  outside  the  ISA  will  be  zero  (i.e,,  the  amount  of  FL/I 
storage  required  and  the  ISA  will  be  equal).  To  establish  the 
correct  value  of  the  ISA,  an  estimate  of  the  proper  size  is  made 
and  a report  is  asked  for  via  the  PARM  field  of  the  EXEC  card  for 
the  job.  Upon  receipt  of  the  report,  the  correct  amount  of  PL/I 
storage  required  will  be  known,  and  subsequent  runs  should  use 
this  value  as  the  ISASIZE.  If  program  variables  that  affect  the 
dynamic  storage  are  changed  between  runs,  the  amount  of  proper 
ISA  storage  may  change.  For  this  reason,  it  is  generally  good 
practice  to  have  storage  reports  made  with  each  run  and  make  ad- 
justments to  the  ISASIZE  as  required. 


INTERFACE  WITH  THE  OPERATING  SYSTEM 

The  machine  configuration  for  the  URLIM  program  code  and 
the  SHTP  subroutines  is  the  IBM  360.  The  PL/I  source  was  compiled 
with  the  optimizer  version  of  the  PL/I  compiler.  All  of  the  rou- 
tines used  by  URLIM  or  contained  in  the  SHTP  library  are  available 
as  load  modules.  Also,  the  URLIM  main  program  is  available  ts  a 
load  module.  The  extent  of  interaction  with  the  IBM  operating 
system  (OS)  is  then  to  properly  assemble  the  required  subroutine 
modules  and  execute  the  program.  The  discussions  that  follow  will 
be  applicable  to  the  running  of  URLIM  or  to  the  execution  of  an- 
other application  program  using  the  SHTP  library,  with  the  follow- 
ing difference.  For  URLIM  runs,  the  main  program  is  already 
written,  compiled,  and  stored  in  a data  set;  for  a run  with  SHTP 
modules,  the  main  program  will  be  user-supplied  and  assumed  avail- 
able as  a precompiled  load  module  ready  for  use. 

To  further  preface  the  following  discussions,  a general 
description  of  the  IBM  360  system  environment  will  be  made.  OS 
is  fundamentally  a supervisory  program  that  oversees  the  illoca- 
tion  of  the  360's  basic  resources,  namely  central  processor  unit 
(CPU)  time,  main  storage  (region),  and  peiipneral  storage  devices. 
Interaction  with  OS  is  accomplished  by  writing  statements  in  job 
control  language  (JCL) . These  statements  serve  as  preparatory 
remarks  to  OS  in  that  they  identify  the  job  to  be  run  and  make 
requests  for  some  of  the  three  aforementioned  basic  resources. 
There  are  three  fundamental  types  of  JCL  statements:  (1)  the  JOB 

card  (containing  user  identification  and  accounting  information) 
that  serves  to  identify  distinct  jobs  to  OS:  (2)  EXEC  cards  that 
identify  what  program  is  to  be  run,  and  (3)  data  definition  (DD) 
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cards  that  identify  the  various  files  (data-sets)  that  will  be  re- 
quired by  the  program.  Requests  for  region  (addressable  storage) 
and  CPU  time  are  made  on  the  EXEC  card. 

At  any  given  360  installation  there  will  be  a set  of  pro- 
grams that  are  used  very  frequently  by  many  users.  The  attendant 
JCL  (which  may  be  quite  lengthy)  is  often  contained  in  files  known 
as  "cataloged  procedures."  These  high-usage  programs  generally 
include  the  IBM-supplied  high-level  compilers  (FL/I,  Fortran,  etc.), 
the  link-editor  program,  and  others.  The  JCL  described  here  will 
include  (1)  that  which  is  needed  to  use  the  link-editor  and 
(2)  that  which  is  sufficient  to  run  a previously  link-edited  module. 

The  Link-Editor 


The  cataloged  procedure  OL  used  at  the  IBM  360/91  installa- 
tion at  APL  is  used  to  execute  the  link-editor;  OLG  is  used  to 
execute  the  link-editor  and  then  execute  the  resultant  program. 

The  use  of  procedure  OL  is  shown  in  the  example  below.  OL  is  used 
to  assemble  the  required  load  modules  together  and  save  the  re- 
sultant, fully  link-edited  program  on  the  file  described  by  the 
L.SYSLMOD  DD  card: 

//  EXEC  OL,'LIB-BBE. FRAZER. BASIC02', 

//  PARM . L- ' LI ST ,MAP , LET ' 

//L.SYSLMOD  DD  DSN-xxx. xxx.xxx(yyy) , 

//  DISP-(NEW,CATLG), 

//  SPACE- (3156, (60, 20, 2), RLSE), 

//  UNIT-SAVE 

//L.SYSIN  DD  * 

INCLUDE  SYSLIB  (URLIM) 

ENTRY  PLISTART 

/* 

The  example  first  names  the  cataloged  data  set  BBK. FRAZER. 
BASIC02  as  the  SYSLIB  file  for  the  link-editor.  Next,  the  file 
SYSLMOD  is  designated  to  be  data  set  xxx.xxx.xxx(yyy)  (i.e., 
member  yyy  of  the  partitioned  data  set  xxx.xxx.xxx) . Further,  the 
data  set  xxx.xxx.xxx  is  to  be  saved  via  cataloging  for  later  use. 
The  "PARM.L-"  specifies  operating  conditions  for  the  link-editor, 
namely,  that  a storage  map  is  requested  (MAP)  and  will  be  printed 
by  the  link-editor;  that  all  input  to  the  link-editor  will  be 
listed  (LIST);  that  the  link-editor  shall  continue  its  operation 
even  though  an  error  may  occur  (LET). 

The  input  file  to  the  link-editor  program  (//L.SYSIN)  con- 
tains two  instructions  to  the  link-editor:  (1)  to  include  the  mem- 
ber URLIM  from  the  SYSLIB  data  set  and  (2)  that  the  entrv  point 
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(first  executable  statement)  of  the  module  is  named  PLISTART,  a 
name  provided  automatically  by  the  PL/I  optimizer  compiler  and 
found  with  the  main  program.  The  member  URLIM  in  this  example  is 
the  main  program  and  will  be  looked  for  by  the  link  editor  as  a 
member  of  the  partitioned  data  set  described  by  the  SYSLIB  DD  card 
in  the  cataloged  procedure  OL.  All  subroutines  called  by  the 
URLIM  program  and  all  subroutines  subsequently  called  are  assumed 
to  be  members  of  either  the  data  set  BBE. FRAZER. BAS IC02,  or  of 
the  system  libraries  named  automatically  within  the  cataloged  pro- 
cedure OL.  The  link-editor's  output  (the  complete  executable  pro- 
gram) is  written  on  the  file  SYSLMOD  and  the  job  is  finished.  The 
SYSLMOD  definitions  indicate  to  the  system  that  the  file  is  to  be 
saved  and  cata/.ged  for  later  use. 

The  interested  reader  will  find  further  descriptions  of  the 
link-editor  program  and  descriptions  of  other  ways  of  assembling 
a complete  executable  program  in  Refs.  16  and  17.  These  other 
methods  include  overlay  defining  to  maximize  storage  use  and  ex- 
cluding routines  that  are  referred  to  but  never  actually  called. 


The  procedure  OLG  is  used  in  the  same  way  as  procedure  OL 
except  that  the  completely  executable  program  (stored  on  the 
SYSLMOD  file)  is  executed  in  a subsequent  JOB  step: 


OLG, LIB-'BBE. FRAZER. BASIC02' , 
PARM.L-’ LIST, MAP, LET', 
PARM.G-'ISASIZE(iiK),R' 

DD  DSN" xx>, . xxx . xxx  (yyy ) , 
DISP-(NEW.CATLG) , UNIT-SAVE, 
SPACE- (3156,(60,20,2), RLSE) 

A 

SYSLIB(URLIM) 

PLISTART 


//  EXEC 

// 

// 

//L. SYSLMOD 

// 

// 

//L.SYSIN  DD 
INCLUDE 
ENTRY 

/* 

//G.SYSIN  DD  * 

(input  data,  as  required) 
//G.READFIL  DD  DSN-xyz .abc .DISP-SHR 


The  PARM.G  statement  passes  parameters  to  the  executing  pro- 
gram and  are  described  in  detail  in  the  preceding  discussion.  The 
input  files  required  by  the  program  are  included  at  the  and  of  the 


Ref,  16.  ‘'System  360  Model  9i  User's  Guide,"  APL/JHU  BCS-1:A0, 
November  1973. 

Ref.  17.  "IBM  OS  Linkage  Editor  and  Loader,"  Eleventh  Edition, 
IBM  File  No.  S360/S370-31,  Order  No.  GC28-6538-10,  April  1973. 
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JCL  prefixed  with  "G.:'  to  ir.dicate  the  second  job  step  (named 
"G")  in  the  same  way  as  the  xink-editor  step  files  are  prefixed 
with  "L." 


Complete  Executable  Programs 


Once  a complete  executable  program  has  been  created  (as  de" 
scribed  above  or  by  the  use  of  other  cataloged  procedures) , an 
execution  may  be  obtained  by  providing  OS  with  the  proper  time  and 
region  requirements  as  well  as  the  required  file  definitions  (via 
DD  cards) . The  JCL  given  below  is  an  example  of  what  can  be  spe- 
cified: 


//  EXEC  PGM-URLIM, 

//  REGION-rrrK, 

//  TIME-tt, 

//  PARM-' ISASIZE(iiiK) , REPORT' 

//STEPLIB  Di)  DSN-xxx.xxx.xxx,DISP-SHR 

//PLIDUMP  DD  SYSOUT-=A 

//SYSIN  DD  * 


(input  data  to  follow  here) 

/* 

//READFIL  DD  DSN-xyt.abc.DISP-SHR 


In  this  example,  the  EXEC  card  specifies  (via  the  PGM-  state 
ment)  the  member  within  the  PDS  viamed  on  the  STEPLIB  card  (in  this 
case  the  data  set  named  xxx.xxx.xxx)  that  contains  the  completely 
executable  program.  (This  example  Is  consistent  with  the  examples 
above  for  OL  and  OLG  in  that  the  names  coincide.)  The  EXEC  card 
asks  for  rrrK  (rrr-thousand)  bytes  ("characters")  of  storage  and 
requests  tt  minutes  of  CPU  time.  [The  time  request  can  be  speci- 
fied in  minutes  and  seconds  as  "TIME- (mm, ss)"] . Additionally,  a 
parameter  is  passed  to  the  executing  program  as  "PARM-...."  This 
parameter  specifies  that  the  Initial  Storage  Area  (ISA)  is  to  be 
iiiK  bytes  long  and  that  a storage  report  is  to  be  given.  The 
significance  of  the  ISASIZE,  and  the  interpretation  of  the  REPORT 
are  discussed  above  in  the  section  on  Dynamic  Storage.  The 
STEPLIB  card  names  the  data  set  containing  the  program  to  be  exe- 
cuted (URLIM  in  this  instance).  The  PLIDUMP  file  is  the  print  file 
onto  which  the  REPORT  will  be  written.  The  SYSIN  file  includes 
inpvt  required  by  the  program  and  expected  from  the  file  SYSIN. 

The  READFIL  file  is  indicative  of  how  files  other  than  SYSIN  can 
be  used  for  supplying  input  data  to  the  various  READ  routines  (e.g. 
READRK) . These  routines,  according  to  values  supplied  as  argu- 
ments, can  read  the  required  data  from  any  file  named  by  a DD 
statement  included  with  the  JCL  fc  the  execution  of  the  job;  the 
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card  "//READFIL  DD  above  describes  such  a file.  Other  spe- 

cific examples  of  these  auxiliary  files  are  given  in  the  READ  rou- 
tine usage  descriptions  m Vol.  2 of  this  report. 

This  demonstrates  the  way  in  which  an  URLIM  or  SHTP  program 
module  can  be  assembled  and  executed  in  an  IBM  OS  360  system  en- 
vironment. There  are  other  ways  of  accomplishing  the  same  results 
and  the  experienced  user  will  experiment  with  and  exploit  such 
avenues  that  prove  beneficial  to  his  needs. 
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LIST  OF  SYMBOLS 

Symbol  Definition  and  Assumed  Units 

A area  (ft^) 

a sonic  velocity  (ft/s) 

radial  position 

C vector  of  heat  capacity  values  for  each  node  (Btu/ft^-BR) 

C , C normal  and  axial  force  coefficients  for  the  radome 
Nr  Ar 

C^,  C2  constants  in  Planck's  equation  (Eq . (23) ) and  others 
c^  friction  coefficient 

cv  specific  heat  (Btu/lbm-°R) , equal  to  du/dt 

0^  specific  heat  (Btu/lbm-°R) 

d wall  thickness  or  length  measure  (in.) 

D antenna  diameter  (in,);  arbitrary  constant 

E,  black  body  irradiance  (Btu/ft^-s) 

b 3 

E grey  body  irradiance  (Btu/ft  -s) 

8 2 
E^  Young '8  Modulus  of  the  ith  subregion  (lb/in  ) 

Ae/e^  boresight  error  slope  (deg/deg) 

F^r*  F normal  and  axial  forces  acting  on  the  radome  shape  ( lbf ) 
Af/fg  percent  change  in  frequency 

F , f resultant  forces  acting  on  radome  in  x and  y 

x ^ directions  (ioi) 

2 

h^  heat  transfer  coefficient  (lbm/ft  -s) 

H the  vector  of  total  conductance  for  each  node  n 

(Btu/h-#R) 

I value  of  the  general  integral  J F(x)dx 

i enthalpy  (Btu/lbm) 

k thermal  conductivity  (Btu/ft~°R-h) 

K the  net  thermal  conductivity  between  two  dissimilar 

tC  materials,  including  contact  resistance  (Btu/ft-°R-h) 

K number  of  nodes  connected  to  node  n via  internal 

conduction 
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L 

M 

m 

r 

M0 

Nu 

Pr 

P,  P 

Q 

q 

• a 

q 

M 

q 

Q* 

R 

R 

R 

*b 

r 

r 

r 

St 


t 

t 

u 

V 

V 
v 

X,  y,  z 
X,  Y 


Z 


length  measure  (ft) 

Mach  number 
radome  mass  (lbm) 

resultant  bending  moment  acting  at  radome's  base  (in.-lbf) 
Nusselt  number,  equal  to  c^h^x/k 
Prandtl  number  equal  to  Pc  /k 

2 P 

pressure  (lbf/ft  ) 

heat  flow  vector  (Btu/h) 

magnitude  of  heat  flow  vector  ( | Q J ) (Btu/h) 

generalized  heat  flux  term  (Btu/ft^-h) 

generalized  heat  generation  term  (Btu/h-ft^) 

the  number  of  independent  heat  flux  sources  exchanging 
heat  with  node  n 

number  of  nodes  connected  to  node  n via  radiation,  or 
relaxation  factor 

gas  constant  (ft-lbf/lbm-®R) 

resultant  force  (lbf) 

radome  base  radius  (in.) 

radial  coordinate  (ft) 

recovery  factor 

reflectivity 

Stanton  number  • q/pV(i  ~ i ) 

r w 

temperature  (*R) 

time  (8  or  h) 

thickness  (in.) 

internal  energy  (Btu/lbm) 

velocity  (ft/s) 

3 

volume  (ft  ) 

3 

specific  volume  (lbra/ft  ) 
spatial  coordinates  (ft) 

coordinate  axes  or  generalized  functional  values 
(i.e.,  X - f (Y) ) 

compressibility  of  a gas 
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i 

I 

I 

I 

I 

I 

1 


Greek  Symbols 


a 

a 

a,  a , 8,  8 

Y 

V 
6 
e 
e 

0,  ♦ 

6 

X 

u 

v 

a 

T 

a 

z 

°e 

a 

r 

•)> 


absorptivity 

thermal  expansion  coefficient  (°R-^) 

constants  in  heat  transfer  correlations 

ratio  of  specific  heats 

gradient  vector  operator  (ft“*) 

error  correction  term  (°R) 

emissivity 

dielectric  constart 

angular  coordinates  (rad) 

quadrant  elevation  angle  (rad) 

wavelength  (in.) 

fluid  viscosity  (lbm/ft-s) 

Poisson' 8 ratio 

“*8  2 a 

Stephan-Boltzmann  constant  (0.174  * 10  Btu/ft  -h-8R; 

fluid/wall  shear  stress  (lbf/in^) 

2 

axial  direction  stress  (lb/in  ) 

2 

circumferential  direction  stress  (lb/in  ) 

2 

radial  direction  stress  (lb/in  ) 
electrical  thickness  of  radome  wall 


& 


•/it 


■iis 


' &• 
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ERRATA  SHEET  FOR  APL/JHU  TG  1293A, 

URLIM  — A Unified  Radome  Limitations  Computer  Program, 
Volume  1 — Theoretical  Background 


1.  Change  Eq.  35  (p.  30)  to: 


i+l 


i+l 


i+l 


i » 1 to  (n  - 1)  , 


(35) 


2.  Change  the  third  line  of  the  paragraph  following  Eq.  40 
(p.  31)  to: 

stitution  of  these  into  Eqs.  (29)  through  (31)  yields  o , a , and 

^ j ” j 


3.  Change  the  first  line  of  the  paragraph  following  Eq.  51c 
(p.  45)  to: 

The  stresses  caused  by  these  forces  are  distributed  around 

4.  Change  the  equation  on  the  sixth  line  of  p.  61  to: 

-2C/r3 4 5 6  + (2/r)(C/r2)  - 0 , 

5.  Change  the  third  line  following  Eq.  69  (p.  66)  to: 

T.  Combining  Eqs.  (68)  and  (69)  yields: 

6.  Change  the  third  line  of  the  last  paragraph  on  p.  66  to: 
correction  factors  (Qq/Q'q)*  The  heating  rate  for  each  node  (Q) 
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7.  Change  Eq.  75  (p.  68)  to: 


3k 


n-1 


3T 


f'l  [(Tn  + Ti)/23 
2 


(75) 


8.  Change  the  fourth  line  of  the  second  paragraph  following 
Eq.  76  (p.  68)  to: 

on  the  initial  guess  for  the  temperature  field  (TQ);  (2)  the 

9.  Change  the  first  line  of  the  second  paragraph  on  p.  70  to: 

The  above  system  shows  nodes  1 through  n as  regular  inter- 

10.  Change  the  eighth  line  on  p.  74  to: 

is  the  extent  of  the  vector  then  i is  set  to  n-1.  Exactly 

11.  Change  the  fifteenth  line  of  the  second  paragraph  on  p.  79  to: 
age  (c.f.  Appendix  F of  Vol.  2 of  this  report).  When  employing 

12.  Change  the  last  line  of  the  second  paragraph  on  p.  79  to: 

Vol.  2 of  this  report). 

13.  Change  the  seventh  line  on  p.  87  to: 

R.  W.  Newman  and  D.  Brockelbank  for  their  suggestions  and  assistance 


